We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
griddata linear / LinearNDInterpolator unexpected behavior
I was trying to speed up my code recently, by reducing the data I feed into griddata with method linear or in other words LinearNDInterpolator. I was greatly surprised to find out that this altered my interpolated values significantly.

Reading the docs of LinearNDInterpolator my understanding of this LinearNDInterpolator or linear method of griddata respectivly is the following:
The surface is modelled with triangular planes between the points. Interpolation to a certain new coordinate should be only dependent on the three next nearest coordinates and values.
Therefore, in principle I should be able to only keep the three nearest data points for each coordinate of interest. Especially it shouldn't matter if I remove data that is not anywhere near the area of concern.

I wrote a "minimal example" that illustrates the problem well. 
1. I define an outer grid with 10x10 positions and random values between 0 and 1.
2. From that I copy a smaller inner grid with the exact same values with size 5x5.
3. I define 50 random positions within this inner grid.
4. I perform liner interpolation to the coordinates defined in 3. based on maps defined in 1. and 2.
5. I measure the difference between the two interpolations.

For many points the two interpolations are identical. However, there is always a difference for some positions with significant value like .2 or .5.

```
import numpy as np
import pandas as pd
from scipy.interpolate import griddata

low, high = -10, 10
low_in, high_in = -2, 2

points_huge = pd.DataFrame(
    [[i, j, np.random.rand()] for i in range(low, high) for j in range(low, high)],
    columns=['x', 'y', 'value']
)

points_in = points_huge[
    (points_huge.x >= low_in) &
    (points_huge.x <= high_in) &
    (points_huge.y >= low_in) &
    (points_huge.y <= high_in)
].copy()

eval_rand = pd.DataFrame(
    [[
        low_in + (high_in - low_in) * np.random.rand() * 0.9,
        low_in + (high_in - low_in) * np.random.rand() * 0.9
    ] for i in range(50)],
    columns=['x', 'y']
)

eval_rand['value_huge'] = griddata(
    points_huge[['x', 'y']],
    points_huge.value,
    eval_rand[['x', 'y']],
    method='linear'
)

eval_rand['value_in'] = griddata(
    points_in[['x', 'y']],
    points_in.value,
    eval_rand[['x', 'y']],
    method='linear'
)

print(max(abs(eval_rand.value_huge-eval_rand.value_in)))
```
I saw that the inner interpolation is done by QHull and found on the page of QHull that triangulation of non-convex surfaces is not supported. Does that mean my random surface can't be interpolated correctly?
![Capture](https://user-images.githubusercontent.com/26696299/200874003-ae8b95bc-77d8-465e-bec3-676a64dc304f.PNG)


--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `scipy/interpolate/_ndgriddata.py`
Content:
```
1 """
2 Convenience interface to N-D interpolation
3 
4 .. versionadded:: 0.9
5 
6 """
7 import numpy as np
8 from .interpnd import LinearNDInterpolator, NDInterpolatorBase, \
9      CloughTocher2DInterpolator, _ndim_coords_from_arrays
10 from scipy.spatial import cKDTree
11 
12 __all__ = ['griddata', 'NearestNDInterpolator', 'LinearNDInterpolator',
13            'CloughTocher2DInterpolator']
14 
15 #------------------------------------------------------------------------------
16 # Nearest-neighbor interpolation
17 #------------------------------------------------------------------------------
18 
19 
20 class NearestNDInterpolator(NDInterpolatorBase):
21     """NearestNDInterpolator(x, y).
22 
23     Nearest-neighbor interpolation in N > 1 dimensions.
24 
25     .. versionadded:: 0.9
26 
27     Methods
28     -------
29     __call__
30 
31     Parameters
32     ----------
33     x : (Npoints, Ndims) ndarray of floats
34         Data point coordinates.
35     y : (Npoints,) ndarray of float or complex
36         Data values.
37     rescale : boolean, optional
38         Rescale points to unit cube before performing interpolation.
39         This is useful if some of the input dimensions have
40         incommensurable units and differ by many orders of magnitude.
41 
42         .. versionadded:: 0.14.0
43     tree_options : dict, optional
44         Options passed to the underlying ``cKDTree``.
45 
46         .. versionadded:: 0.17.0
47 
48 
49     Notes
50     -----
51     Uses ``scipy.spatial.cKDTree``
52 
53     Examples
54     --------
55     We can interpolate values on a 2D plane:
56 
57     >>> from scipy.interpolate import NearestNDInterpolator
58     >>> import numpy as np
59     >>> import matplotlib.pyplot as plt
60     >>> rng = np.random.default_rng()
61     >>> x = rng.random(10) - 0.5
62     >>> y = rng.random(10) - 0.5
63     >>> z = np.hypot(x, y)
64     >>> X = np.linspace(min(x), max(x))
65     >>> Y = np.linspace(min(y), max(y))
66     >>> X, Y = np.meshgrid(X, Y)  # 2D grid for interpolation
67     >>> interp = NearestNDInterpolator(list(zip(x, y)), z)
68     >>> Z = interp(X, Y)
69     >>> plt.pcolormesh(X, Y, Z, shading='auto')
70     >>> plt.plot(x, y, "ok", label="input point")
71     >>> plt.legend()
72     >>> plt.colorbar()
73     >>> plt.axis("equal")
74     >>> plt.show()
75 
76     See also
77     --------
78     griddata :
79         Interpolate unstructured D-D data.
80     LinearNDInterpolator :
81         Piecewise linear interpolant in N dimensions.
82     CloughTocher2DInterpolator :
83         Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D.
84 
85     """
86 
87     def __init__(self, x, y, rescale=False, tree_options=None):
88         NDInterpolatorBase.__init__(self, x, y, rescale=rescale,
89                                     need_contiguous=False,
90                                     need_values=False)
91         if tree_options is None:
92             tree_options = dict()
93         self.tree = cKDTree(self.points, **tree_options)
94         self.values = np.asarray(y)
95 
96     def __call__(self, *args):
97         """
98         Evaluate interpolator at given points.
99 
100         Parameters
101         ----------
102         x1, x2, ... xn : array-like of float
103             Points where to interpolate data at.
104             x1, x2, ... xn can be array-like of float with broadcastable shape.
105             or x1 can be array-like of float with shape ``(..., ndim)``
106 
107         """
108         xi = _ndim_coords_from_arrays(args, ndim=self.points.shape[1])
109         xi = self._check_call_shape(xi)
110         xi = self._scale_x(xi)
111         dist, i = self.tree.query(xi)
112         return self.values[i]
113 
114 
115 #------------------------------------------------------------------------------
116 # Convenience interface function
117 #------------------------------------------------------------------------------
118 
119 def griddata(points, values, xi, method='linear', fill_value=np.nan,
120              rescale=False):
121     """
122     Interpolate unstructured D-D data.
123 
124     Parameters
125     ----------
126     points : 2-D ndarray of floats with shape (n, D), or length D tuple of 1-D ndarrays with shape (n,).
127         Data point coordinates.
128     values : ndarray of float or complex, shape (n,)
129         Data values.
130     xi : 2-D ndarray of floats with shape (m, D), or length D tuple of ndarrays broadcastable to the same shape.
131         Points at which to interpolate data.
132     method : {'linear', 'nearest', 'cubic'}, optional
133         Method of interpolation. One of
134 
135         ``nearest``
136           return the value at the data point closest to
137           the point of interpolation. See `NearestNDInterpolator` for
138           more details.
139 
140         ``linear``
141           tessellate the input point set to N-D
142           simplices, and interpolate linearly on each simplex. See
143           `LinearNDInterpolator` for more details.
144 
145         ``cubic`` (1-D)
146           return the value determined from a cubic
147           spline.
148 
149         ``cubic`` (2-D)
150           return the value determined from a
151           piecewise cubic, continuously differentiable (C1), and
152           approximately curvature-minimizing polynomial surface. See
153           `CloughTocher2DInterpolator` for more details.
154     fill_value : float, optional
155         Value used to fill in for requested points outside of the
156         convex hull of the input points. If not provided, then the
157         default is ``nan``. This option has no effect for the
158         'nearest' method.
159     rescale : bool, optional
160         Rescale points to unit cube before performing interpolation.
161         This is useful if some of the input dimensions have
162         incommensurable units and differ by many orders of magnitude.
163 
164         .. versionadded:: 0.14.0
165 
166     Returns
167     -------
168     ndarray
169         Array of interpolated values.
170 
171     Notes
172     -----
173 
174     .. versionadded:: 0.9
175 
176     For data on a regular grid use `interpn` instead.
177 
178     Examples
179     --------
180 
181     Suppose we want to interpolate the 2-D function
182 
183     >>> import numpy as np
184     >>> def func(x, y):
185     ...     return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2
186 
187     on a grid in [0, 1]x[0, 1]
188 
189     >>> grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]
190 
191     but we only know its values at 1000 data points:
192 
193     >>> rng = np.random.default_rng()
194     >>> points = rng.random((1000, 2))
195     >>> values = func(points[:,0], points[:,1])
196 
197     This can be done with `griddata` -- below we try out all of the
198     interpolation methods:
199 
200     >>> from scipy.interpolate import griddata
201     >>> grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')
202     >>> grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')
203     >>> grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic')
204 
205     One can see that the exact result is reproduced by all of the
206     methods to some degree, but for this smooth function the piecewise
207     cubic interpolant gives the best results:
208 
209     >>> import matplotlib.pyplot as plt
210     >>> plt.subplot(221)
211     >>> plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')
212     >>> plt.plot(points[:,0], points[:,1], 'k.', ms=1)
213     >>> plt.title('Original')
214     >>> plt.subplot(222)
215     >>> plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')
216     >>> plt.title('Nearest')
217     >>> plt.subplot(223)
218     >>> plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')
219     >>> plt.title('Linear')
220     >>> plt.subplot(224)
221     >>> plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')
222     >>> plt.title('Cubic')
223     >>> plt.gcf().set_size_inches(6, 6)
224     >>> plt.show()
225 
226     See Also
227     --------
228     LinearNDInterpolator :
229         Piecewise linear interpolant in N dimensions.
230     NearestNDInterpolator :
231         Nearest-neighbor interpolation in N dimensions.
232     CloughTocher2DInterpolator :
233         Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D.
234 
235     """
236 
237     points = _ndim_coords_from_arrays(points)
238 
239     if points.ndim < 2:
240         ndim = points.ndim
241     else:
242         ndim = points.shape[-1]
243 
244     if ndim == 1 and method in ('nearest', 'linear', 'cubic'):
245         from ._interpolate import interp1d
246         points = points.ravel()
247         if isinstance(xi, tuple):
248             if len(xi) != 1:
249                 raise ValueError("invalid number of dimensions in xi")
250             xi, = xi
251         # Sort points/values together, necessary as input for interp1d
252         idx = np.argsort(points)
253         points = points[idx]
254         values = values[idx]
255         if method == 'nearest':
256             fill_value = 'extrapolate'
257         ip = interp1d(points, values, kind=method, axis=0, bounds_error=False,
258                       fill_value=fill_value)
259         return ip(xi)
260     elif method == 'nearest':
261         ip = NearestNDInterpolator(points, values, rescale=rescale)
262         return ip(xi)
263     elif method == 'linear':
264         ip = LinearNDInterpolator(points, values, fill_value=fill_value,
265                                   rescale=rescale)
266         return ip(xi)
267     elif method == 'cubic' and ndim == 2:
268         ip = CloughTocher2DInterpolator(points, values, fill_value=fill_value,
269                                         rescale=rescale)
270         return ip(xi)
271     else:
272         raise ValueError("Unknown interpolation method %r for "
273                          "%d dimensional data" % (method, ndim))
274
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/scipy/interpolate/_ndgriddata.py b/scipy/interpolate/_ndgriddata.py
--- a/scipy/interpolate/_ndgriddata.py
+++ b/scipy/interpolate/_ndgriddata.py
@@ -50,6 +50,8 @@
     -----
     Uses ``scipy.spatial.cKDTree``
 
+    .. note:: For data on a regular grid use `interpn` instead.
+
     Examples
     --------
     We can interpolate values on a 2D plane:
@@ -173,7 +175,7 @@
 
     .. versionadded:: 0.9
 
-    For data on a regular grid use `interpn` instead.
+    .. note:: For data on a regular grid use `interpn` instead.
 
     Examples
     --------

{"golden_diff": "diff --git a/scipy/interpolate/_ndgriddata.py b/scipy/interpolate/_ndgriddata.py\n--- a/scipy/interpolate/_ndgriddata.py\n+++ b/scipy/interpolate/_ndgriddata.py\n@@ -50,6 +50,8 @@\n     -----\n     Uses ``scipy.spatial.cKDTree``\n \n+    .. note:: For data on a regular grid use `interpn` instead.\n+\n     Examples\n     --------\n     We can interpolate values on a 2D plane:\n@@ -173,7 +175,7 @@\n \n     .. versionadded:: 0.9\n \n-    For data on a regular grid use `interpn` instead.\n+    .. note:: For data on a regular grid use `interpn` instead.\n \n     Examples\n     --------\n", "issue": "griddata linear / LinearNDInterpolator unexpected behavior\nI was trying to speed up my code recently, by reducing the data I feed into griddata with method linear or in other words LinearNDInterpolator. I was greatly surprised to find out that this altered my interpolated values significantly.\r\n\r\nReading the docs of LinearNDInterpolator my understanding of this LinearNDInterpolator or linear method of griddata respectivly is the following:\r\nThe surface is modelled with triangular planes between the points. Interpolation to a certain new coordinate should be only dependent on the three next nearest coordinates and values.\r\nTherefore, in principle I should be able to only keep the three nearest data points for each coordinate of interest. Especially it shouldn't matter if I remove data that is not anywhere near the area of concern.\r\n\r\nI wrote a \"minimal example\" that illustrates the problem well. \r\n1. I define an outer grid with 10x10 positions and random values between 0 and 1.\r\n2. From that I copy a smaller inner grid with the exact same values with size 5x5.\r\n3. I define 50 random positions within this inner grid.\r\n4. I perform liner interpolation to the coordinates defined in 3. based on maps defined in 1. and 2.\r\n5. I measure the difference between the two interpolations.\r\n\r\nFor many points the two interpolations are identical. However, there is always a difference for some positions with significant value like .2 or .5.\r\n\r\n```\r\nimport numpy as np\r\nimport pandas as pd\r\nfrom scipy.interpolate import griddata\r\n\r\nlow, high = -10, 10\r\nlow_in, high_in = -2, 2\r\n\r\npoints_huge = pd.DataFrame(\r\n    [[i, j, np.random.rand()] for i in range(low, high) for j in range(low, high)],\r\n    columns=['x', 'y', 'value']\r\n)\r\n\r\npoints_in = points_huge[\r\n    (points_huge.x >= low_in) &\r\n    (points_huge.x <= high_in) &\r\n    (points_huge.y >= low_in) &\r\n    (points_huge.y <= high_in)\r\n].copy()\r\n\r\neval_rand = pd.DataFrame(\r\n    [[\r\n        low_in + (high_in - low_in) * np.random.rand() * 0.9,\r\n        low_in + (high_in - low_in) * np.random.rand() * 0.9\r\n    ] for i in range(50)],\r\n    columns=['x', 'y']\r\n)\r\n\r\neval_rand['value_huge'] = griddata(\r\n    points_huge[['x', 'y']],\r\n    points_huge.value,\r\n    eval_rand[['x', 'y']],\r\n    method='linear'\r\n)\r\n\r\neval_rand['value_in'] = griddata(\r\n    points_in[['x', 'y']],\r\n    points_in.value,\r\n    eval_rand[['x', 'y']],\r\n    method='linear'\r\n)\r\n\r\nprint(max(abs(eval_rand.value_huge-eval_rand.value_in)))\r\n```\r\nI saw that the inner interpolation is done by QHull and found on the page of QHull that triangulation of non-convex surfaces is not supported. Does that mean my random surface can't be interpolated correctly?\r\n![Capture](https://user-images.githubusercontent.com/26696299/200874003-ae8b95bc-77d8-465e-bec3-676a64dc304f.PNG)\r\n\n", "before_files": [{"content": "\"\"\"\nConvenience interface to N-D interpolation\n\n.. versionadded:: 0.9\n\n\"\"\"\nimport numpy as np\nfrom .interpnd import LinearNDInterpolator, NDInterpolatorBase, \\\n     CloughTocher2DInterpolator, _ndim_coords_from_arrays\nfrom scipy.spatial import cKDTree\n\n__all__ = ['griddata', 'NearestNDInterpolator', 'LinearNDInterpolator',\n           'CloughTocher2DInterpolator']\n\n#------------------------------------------------------------------------------\n# Nearest-neighbor interpolation\n#------------------------------------------------------------------------------\n\n\nclass NearestNDInterpolator(NDInterpolatorBase):\n    \"\"\"NearestNDInterpolator(x, y).\n\n    Nearest-neighbor interpolation in N > 1 dimensions.\n\n    .. versionadded:: 0.9\n\n    Methods\n    -------\n    __call__\n\n    Parameters\n    ----------\n    x : (Npoints, Ndims) ndarray of floats\n        Data point coordinates.\n    y : (Npoints,) ndarray of float or complex\n        Data values.\n    rescale : boolean, optional\n        Rescale points to unit cube before performing interpolation.\n        This is useful if some of the input dimensions have\n        incommensurable units and differ by many orders of magnitude.\n\n        .. versionadded:: 0.14.0\n    tree_options : dict, optional\n        Options passed to the underlying ``cKDTree``.\n\n        .. versionadded:: 0.17.0\n\n\n    Notes\n    -----\n    Uses ``scipy.spatial.cKDTree``\n\n    Examples\n    --------\n    We can interpolate values on a 2D plane:\n\n    >>> from scipy.interpolate import NearestNDInterpolator\n    >>> import numpy as np\n    >>> import matplotlib.pyplot as plt\n    >>> rng = np.random.default_rng()\n    >>> x = rng.random(10) - 0.5\n    >>> y = rng.random(10) - 0.5\n    >>> z = np.hypot(x, y)\n    >>> X = np.linspace(min(x), max(x))\n    >>> Y = np.linspace(min(y), max(y))\n    >>> X, Y = np.meshgrid(X, Y)  # 2D grid for interpolation\n    >>> interp = NearestNDInterpolator(list(zip(x, y)), z)\n    >>> Z = interp(X, Y)\n    >>> plt.pcolormesh(X, Y, Z, shading='auto')\n    >>> plt.plot(x, y, \"ok\", label=\"input point\")\n    >>> plt.legend()\n    >>> plt.colorbar()\n    >>> plt.axis(\"equal\")\n    >>> plt.show()\n\n    See also\n    --------\n    griddata :\n        Interpolate unstructured D-D data.\n    LinearNDInterpolator :\n        Piecewise linear interpolant in N dimensions.\n    CloughTocher2DInterpolator :\n        Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D.\n\n    \"\"\"\n\n    def __init__(self, x, y, rescale=False, tree_options=None):\n        NDInterpolatorBase.__init__(self, x, y, rescale=rescale,\n                                    need_contiguous=False,\n                                    need_values=False)\n        if tree_options is None:\n            tree_options = dict()\n        self.tree = cKDTree(self.points, **tree_options)\n        self.values = np.asarray(y)\n\n    def __call__(self, *args):\n        \"\"\"\n        Evaluate interpolator at given points.\n\n        Parameters\n        ----------\n        x1, x2, ... xn : array-like of float\n            Points where to interpolate data at.\n            x1, x2, ... xn can be array-like of float with broadcastable shape.\n            or x1 can be array-like of float with shape ``(..., ndim)``\n\n        \"\"\"\n        xi = _ndim_coords_from_arrays(args, ndim=self.points.shape[1])\n        xi = self._check_call_shape(xi)\n        xi = self._scale_x(xi)\n        dist, i = self.tree.query(xi)\n        return self.values[i]\n\n\n#------------------------------------------------------------------------------\n# Convenience interface function\n#------------------------------------------------------------------------------\n\ndef griddata(points, values, xi, method='linear', fill_value=np.nan,\n             rescale=False):\n    \"\"\"\n    Interpolate unstructured D-D data.\n\n    Parameters\n    ----------\n    points : 2-D ndarray of floats with shape (n, D), or length D tuple of 1-D ndarrays with shape (n,).\n        Data point coordinates.\n    values : ndarray of float or complex, shape (n,)\n        Data values.\n    xi : 2-D ndarray of floats with shape (m, D), or length D tuple of ndarrays broadcastable to the same shape.\n        Points at which to interpolate data.\n    method : {'linear', 'nearest', 'cubic'}, optional\n        Method of interpolation. One of\n\n        ``nearest``\n          return the value at the data point closest to\n          the point of interpolation. See `NearestNDInterpolator` for\n          more details.\n\n        ``linear``\n          tessellate the input point set to N-D\n          simplices, and interpolate linearly on each simplex. See\n          `LinearNDInterpolator` for more details.\n\n        ``cubic`` (1-D)\n          return the value determined from a cubic\n          spline.\n\n        ``cubic`` (2-D)\n          return the value determined from a\n          piecewise cubic, continuously differentiable (C1), and\n          approximately curvature-minimizing polynomial surface. See\n          `CloughTocher2DInterpolator` for more details.\n    fill_value : float, optional\n        Value used to fill in for requested points outside of the\n        convex hull of the input points. If not provided, then the\n        default is ``nan``. This option has no effect for the\n        'nearest' method.\n    rescale : bool, optional\n        Rescale points to unit cube before performing interpolation.\n        This is useful if some of the input dimensions have\n        incommensurable units and differ by many orders of magnitude.\n\n        .. versionadded:: 0.14.0\n\n    Returns\n    -------\n    ndarray\n        Array of interpolated values.\n\n    Notes\n    -----\n\n    .. versionadded:: 0.9\n\n    For data on a regular grid use `interpn` instead.\n\n    Examples\n    --------\n\n    Suppose we want to interpolate the 2-D function\n\n    >>> import numpy as np\n    >>> def func(x, y):\n    ...     return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2\n\n    on a grid in [0, 1]x[0, 1]\n\n    >>> grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]\n\n    but we only know its values at 1000 data points:\n\n    >>> rng = np.random.default_rng()\n    >>> points = rng.random((1000, 2))\n    >>> values = func(points[:,0], points[:,1])\n\n    This can be done with `griddata` -- below we try out all of the\n    interpolation methods:\n\n    >>> from scipy.interpolate import griddata\n    >>> grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')\n    >>> grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')\n    >>> grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic')\n\n    One can see that the exact result is reproduced by all of the\n    methods to some degree, but for this smooth function the piecewise\n    cubic interpolant gives the best results:\n\n    >>> import matplotlib.pyplot as plt\n    >>> plt.subplot(221)\n    >>> plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')\n    >>> plt.plot(points[:,0], points[:,1], 'k.', ms=1)\n    >>> plt.title('Original')\n    >>> plt.subplot(222)\n    >>> plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')\n    >>> plt.title('Nearest')\n    >>> plt.subplot(223)\n    >>> plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')\n    >>> plt.title('Linear')\n    >>> plt.subplot(224)\n    >>> plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')\n    >>> plt.title('Cubic')\n    >>> plt.gcf().set_size_inches(6, 6)\n    >>> plt.show()\n\n    See Also\n    --------\n    LinearNDInterpolator :\n        Piecewise linear interpolant in N dimensions.\n    NearestNDInterpolator :\n        Nearest-neighbor interpolation in N dimensions.\n    CloughTocher2DInterpolator :\n        Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D.\n\n    \"\"\"\n\n    points = _ndim_coords_from_arrays(points)\n\n    if points.ndim < 2:\n        ndim = points.ndim\n    else:\n        ndim = points.shape[-1]\n\n    if ndim == 1 and method in ('nearest', 'linear', 'cubic'):\n        from ._interpolate import interp1d\n        points = points.ravel()\n        if isinstance(xi, tuple):\n            if len(xi) != 1:\n                raise ValueError(\"invalid number of dimensions in xi\")\n            xi, = xi\n        # Sort points/values together, necessary as input for interp1d\n        idx = np.argsort(points)\n        points = points[idx]\n        values = values[idx]\n        if method == 'nearest':\n            fill_value = 'extrapolate'\n        ip = interp1d(points, values, kind=method, axis=0, bounds_error=False,\n                      fill_value=fill_value)\n        return ip(xi)\n    elif method == 'nearest':\n        ip = NearestNDInterpolator(points, values, rescale=rescale)\n        return ip(xi)\n    elif method == 'linear':\n        ip = LinearNDInterpolator(points, values, fill_value=fill_value,\n                                  rescale=rescale)\n        return ip(xi)\n    elif method == 'cubic' and ndim == 2:\n        ip = CloughTocher2DInterpolator(points, values, fill_value=fill_value,\n                                        rescale=rescale)\n        return ip(xi)\n    else:\n        raise ValueError(\"Unknown interpolation method %r for \"\n                         \"%d dimensional data\" % (method, ndim))\n", "path": "scipy/interpolate/_ndgriddata.py"}], "after_files": [{"content": "\"\"\"\nConvenience interface to N-D interpolation\n\n.. versionadded:: 0.9\n\n\"\"\"\nimport numpy as np\nfrom .interpnd import LinearNDInterpolator, NDInterpolatorBase, \\\n     CloughTocher2DInterpolator, _ndim_coords_from_arrays\nfrom scipy.spatial import cKDTree\n\n__all__ = ['griddata', 'NearestNDInterpolator', 'LinearNDInterpolator',\n           'CloughTocher2DInterpolator']\n\n#------------------------------------------------------------------------------\n# Nearest-neighbor interpolation\n#------------------------------------------------------------------------------\n\n\nclass NearestNDInterpolator(NDInterpolatorBase):\n    \"\"\"NearestNDInterpolator(x, y).\n\n    Nearest-neighbor interpolation in N > 1 dimensions.\n\n    .. versionadded:: 0.9\n\n    Methods\n    -------\n    __call__\n\n    Parameters\n    ----------\n    x : (Npoints, Ndims) ndarray of floats\n        Data point coordinates.\n    y : (Npoints,) ndarray of float or complex\n        Data values.\n    rescale : boolean, optional\n        Rescale points to unit cube before performing interpolation.\n        This is useful if some of the input dimensions have\n        incommensurable units and differ by many orders of magnitude.\n\n        .. versionadded:: 0.14.0\n    tree_options : dict, optional\n        Options passed to the underlying ``cKDTree``.\n\n        .. versionadded:: 0.17.0\n\n\n    Notes\n    -----\n    Uses ``scipy.spatial.cKDTree``\n\n    .. note:: For data on a regular grid use `interpn` instead.\n\n    Examples\n    --------\n    We can interpolate values on a 2D plane:\n\n    >>> from scipy.interpolate import NearestNDInterpolator\n    >>> import numpy as np\n    >>> import matplotlib.pyplot as plt\n    >>> rng = np.random.default_rng()\n    >>> x = rng.random(10) - 0.5\n    >>> y = rng.random(10) - 0.5\n    >>> z = np.hypot(x, y)\n    >>> X = np.linspace(min(x), max(x))\n    >>> Y = np.linspace(min(y), max(y))\n    >>> X, Y = np.meshgrid(X, Y)  # 2D grid for interpolation\n    >>> interp = NearestNDInterpolator(list(zip(x, y)), z)\n    >>> Z = interp(X, Y)\n    >>> plt.pcolormesh(X, Y, Z, shading='auto')\n    >>> plt.plot(x, y, \"ok\", label=\"input point\")\n    >>> plt.legend()\n    >>> plt.colorbar()\n    >>> plt.axis(\"equal\")\n    >>> plt.show()\n\n    See also\n    --------\n    griddata :\n        Interpolate unstructured D-D data.\n    LinearNDInterpolator :\n        Piecewise linear interpolant in N dimensions.\n    CloughTocher2DInterpolator :\n        Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D.\n\n    \"\"\"\n\n    def __init__(self, x, y, rescale=False, tree_options=None):\n        NDInterpolatorBase.__init__(self, x, y, rescale=rescale,\n                                    need_contiguous=False,\n                                    need_values=False)\n        if tree_options is None:\n            tree_options = dict()\n        self.tree = cKDTree(self.points, **tree_options)\n        self.values = np.asarray(y)\n\n    def __call__(self, *args):\n        \"\"\"\n        Evaluate interpolator at given points.\n\n        Parameters\n        ----------\n        x1, x2, ... xn : array-like of float\n            Points where to interpolate data at.\n            x1, x2, ... xn can be array-like of float with broadcastable shape.\n            or x1 can be array-like of float with shape ``(..., ndim)``\n\n        \"\"\"\n        xi = _ndim_coords_from_arrays(args, ndim=self.points.shape[1])\n        xi = self._check_call_shape(xi)\n        xi = self._scale_x(xi)\n        dist, i = self.tree.query(xi)\n        return self.values[i]\n\n\n#------------------------------------------------------------------------------\n# Convenience interface function\n#------------------------------------------------------------------------------\n\ndef griddata(points, values, xi, method='linear', fill_value=np.nan,\n             rescale=False):\n    \"\"\"\n    Interpolate unstructured D-D data.\n\n    Parameters\n    ----------\n    points : 2-D ndarray of floats with shape (n, D), or length D tuple of 1-D ndarrays with shape (n,).\n        Data point coordinates.\n    values : ndarray of float or complex, shape (n,)\n        Data values.\n    xi : 2-D ndarray of floats with shape (m, D), or length D tuple of ndarrays broadcastable to the same shape.\n        Points at which to interpolate data.\n    method : {'linear', 'nearest', 'cubic'}, optional\n        Method of interpolation. One of\n\n        ``nearest``\n          return the value at the data point closest to\n          the point of interpolation. See `NearestNDInterpolator` for\n          more details.\n\n        ``linear``\n          tessellate the input point set to N-D\n          simplices, and interpolate linearly on each simplex. See\n          `LinearNDInterpolator` for more details.\n\n        ``cubic`` (1-D)\n          return the value determined from a cubic\n          spline.\n\n        ``cubic`` (2-D)\n          return the value determined from a\n          piecewise cubic, continuously differentiable (C1), and\n          approximately curvature-minimizing polynomial surface. See\n          `CloughTocher2DInterpolator` for more details.\n    fill_value : float, optional\n        Value used to fill in for requested points outside of the\n        convex hull of the input points. If not provided, then the\n        default is ``nan``. This option has no effect for the\n        'nearest' method.\n    rescale : bool, optional\n        Rescale points to unit cube before performing interpolation.\n        This is useful if some of the input dimensions have\n        incommensurable units and differ by many orders of magnitude.\n\n        .. versionadded:: 0.14.0\n\n    Returns\n    -------\n    ndarray\n        Array of interpolated values.\n\n    Notes\n    -----\n\n    .. versionadded:: 0.9\n\n    .. note:: For data on a regular grid use `interpn` instead.\n\n    Examples\n    --------\n\n    Suppose we want to interpolate the 2-D function\n\n    >>> import numpy as np\n    >>> def func(x, y):\n    ...     return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2\n\n    on a grid in [0, 1]x[0, 1]\n\n    >>> grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]\n\n    but we only know its values at 1000 data points:\n\n    >>> rng = np.random.default_rng()\n    >>> points = rng.random((1000, 2))\n    >>> values = func(points[:,0], points[:,1])\n\n    This can be done with `griddata` -- below we try out all of the\n    interpolation methods:\n\n    >>> from scipy.interpolate import griddata\n    >>> grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest')\n    >>> grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear')\n    >>> grid_z2 = griddata(points, values, (grid_x, grid_y), method='cubic')\n\n    One can see that the exact result is reproduced by all of the\n    methods to some degree, but for this smooth function the piecewise\n    cubic interpolant gives the best results:\n\n    >>> import matplotlib.pyplot as plt\n    >>> plt.subplot(221)\n    >>> plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')\n    >>> plt.plot(points[:,0], points[:,1], 'k.', ms=1)\n    >>> plt.title('Original')\n    >>> plt.subplot(222)\n    >>> plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')\n    >>> plt.title('Nearest')\n    >>> plt.subplot(223)\n    >>> plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')\n    >>> plt.title('Linear')\n    >>> plt.subplot(224)\n    >>> plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')\n    >>> plt.title('Cubic')\n    >>> plt.gcf().set_size_inches(6, 6)\n    >>> plt.show()\n\n    See Also\n    --------\n    LinearNDInterpolator :\n        Piecewise linear interpolant in N dimensions.\n    NearestNDInterpolator :\n        Nearest-neighbor interpolation in N dimensions.\n    CloughTocher2DInterpolator :\n        Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D.\n\n    \"\"\"\n\n    points = _ndim_coords_from_arrays(points)\n\n    if points.ndim < 2:\n        ndim = points.ndim\n    else:\n        ndim = points.shape[-1]\n\n    if ndim == 1 and method in ('nearest', 'linear', 'cubic'):\n        from ._interpolate import interp1d\n        points = points.ravel()\n        if isinstance(xi, tuple):\n            if len(xi) != 1:\n                raise ValueError(\"invalid number of dimensions in xi\")\n            xi, = xi\n        # Sort points/values together, necessary as input for interp1d\n        idx = np.argsort(points)\n        points = points[idx]\n        values = values[idx]\n        if method == 'nearest':\n            fill_value = 'extrapolate'\n        ip = interp1d(points, values, kind=method, axis=0, bounds_error=False,\n                      fill_value=fill_value)\n        return ip(xi)\n    elif method == 'nearest':\n        ip = NearestNDInterpolator(points, values, rescale=rescale)\n        return ip(xi)\n    elif method == 'linear':\n        ip = LinearNDInterpolator(points, values, fill_value=fill_value,\n                                  rescale=rescale)\n        return ip(xi)\n    elif method == 'cubic' and ndim == 2:\n        ip = CloughTocher2DInterpolator(points, values, fill_value=fill_value,\n                                        rescale=rescale)\n        return ip(xi)\n    else:\n        raise ValueError(\"Unknown interpolation method %r for \"\n                         \"%d dimensional data\" % (method, ndim))\n", "path": "scipy/interpolate/_ndgriddata.py"}]}