test_ndgriddata.py 9.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246
  1. import numpy as np
  2. from numpy.testing import assert_equal, assert_array_equal, assert_allclose
  3. import pytest
  4. from pytest import raises as assert_raises
  5. from scipy.interpolate import (griddata, NearestNDInterpolator,
  6. LinearNDInterpolator,
  7. CloughTocher2DInterpolator)
  8. parametrize_interpolators = pytest.mark.parametrize(
  9. "interpolator", [NearestNDInterpolator, LinearNDInterpolator,
  10. CloughTocher2DInterpolator]
  11. )
  12. class TestGriddata:
  13. def test_fill_value(self):
  14. x = [(0,0), (0,1), (1,0)]
  15. y = [1, 2, 3]
  16. yi = griddata(x, y, [(1,1), (1,2), (0,0)], fill_value=-1)
  17. assert_array_equal(yi, [-1., -1, 1])
  18. yi = griddata(x, y, [(1,1), (1,2), (0,0)])
  19. assert_array_equal(yi, [np.nan, np.nan, 1])
  20. def test_alternative_call(self):
  21. x = np.array([(0,0), (-0.5,-0.5), (-0.5,0.5), (0.5, 0.5), (0.25, 0.3)],
  22. dtype=np.double)
  23. y = (np.arange(x.shape[0], dtype=np.double)[:,None]
  24. + np.array([0,1])[None,:])
  25. for method in ('nearest', 'linear', 'cubic'):
  26. for rescale in (True, False):
  27. msg = repr((method, rescale))
  28. yi = griddata((x[:,0], x[:,1]), y, (x[:,0], x[:,1]), method=method,
  29. rescale=rescale)
  30. assert_allclose(y, yi, atol=1e-14, err_msg=msg)
  31. def test_multivalue_2d(self):
  32. x = np.array([(0,0), (-0.5,-0.5), (-0.5,0.5), (0.5, 0.5), (0.25, 0.3)],
  33. dtype=np.double)
  34. y = (np.arange(x.shape[0], dtype=np.double)[:,None]
  35. + np.array([0,1])[None,:])
  36. for method in ('nearest', 'linear', 'cubic'):
  37. for rescale in (True, False):
  38. msg = repr((method, rescale))
  39. yi = griddata(x, y, x, method=method, rescale=rescale)
  40. assert_allclose(y, yi, atol=1e-14, err_msg=msg)
  41. def test_multipoint_2d(self):
  42. x = np.array([(0,0), (-0.5,-0.5), (-0.5,0.5), (0.5, 0.5), (0.25, 0.3)],
  43. dtype=np.double)
  44. y = np.arange(x.shape[0], dtype=np.double)
  45. xi = x[:,None,:] + np.array([0,0,0])[None,:,None]
  46. for method in ('nearest', 'linear', 'cubic'):
  47. for rescale in (True, False):
  48. msg = repr((method, rescale))
  49. yi = griddata(x, y, xi, method=method, rescale=rescale)
  50. assert_equal(yi.shape, (5, 3), err_msg=msg)
  51. assert_allclose(yi, np.tile(y[:,None], (1, 3)),
  52. atol=1e-14, err_msg=msg)
  53. def test_complex_2d(self):
  54. x = np.array([(0,0), (-0.5,-0.5), (-0.5,0.5), (0.5, 0.5), (0.25, 0.3)],
  55. dtype=np.double)
  56. y = np.arange(x.shape[0], dtype=np.double)
  57. y = y - 2j*y[::-1]
  58. xi = x[:,None,:] + np.array([0,0,0])[None,:,None]
  59. for method in ('nearest', 'linear', 'cubic'):
  60. for rescale in (True, False):
  61. msg = repr((method, rescale))
  62. yi = griddata(x, y, xi, method=method, rescale=rescale)
  63. assert_equal(yi.shape, (5, 3), err_msg=msg)
  64. assert_allclose(yi, np.tile(y[:,None], (1, 3)),
  65. atol=1e-14, err_msg=msg)
  66. def test_1d(self):
  67. x = np.array([1, 2.5, 3, 4.5, 5, 6])
  68. y = np.array([1, 2, 0, 3.9, 2, 1])
  69. for method in ('nearest', 'linear', 'cubic'):
  70. assert_allclose(griddata(x, y, x, method=method), y,
  71. err_msg=method, atol=1e-14)
  72. assert_allclose(griddata(x.reshape(6, 1), y, x, method=method), y,
  73. err_msg=method, atol=1e-14)
  74. assert_allclose(griddata((x,), y, (x,), method=method), y,
  75. err_msg=method, atol=1e-14)
  76. def test_1d_borders(self):
  77. # Test for nearest neighbor case with xi outside
  78. # the range of the values.
  79. x = np.array([1, 2.5, 3, 4.5, 5, 6])
  80. y = np.array([1, 2, 0, 3.9, 2, 1])
  81. xi = np.array([0.9, 6.5])
  82. yi_should = np.array([1.0, 1.0])
  83. method = 'nearest'
  84. assert_allclose(griddata(x, y, xi,
  85. method=method), yi_should,
  86. err_msg=method,
  87. atol=1e-14)
  88. assert_allclose(griddata(x.reshape(6, 1), y, xi,
  89. method=method), yi_should,
  90. err_msg=method,
  91. atol=1e-14)
  92. assert_allclose(griddata((x, ), y, (xi, ),
  93. method=method), yi_should,
  94. err_msg=method,
  95. atol=1e-14)
  96. def test_1d_unsorted(self):
  97. x = np.array([2.5, 1, 4.5, 5, 6, 3])
  98. y = np.array([1, 2, 0, 3.9, 2, 1])
  99. for method in ('nearest', 'linear', 'cubic'):
  100. assert_allclose(griddata(x, y, x, method=method), y,
  101. err_msg=method, atol=1e-10)
  102. assert_allclose(griddata(x.reshape(6, 1), y, x, method=method), y,
  103. err_msg=method, atol=1e-10)
  104. assert_allclose(griddata((x,), y, (x,), method=method), y,
  105. err_msg=method, atol=1e-10)
  106. def test_square_rescale_manual(self):
  107. points = np.array([(0,0), (0,100), (10,100), (10,0), (1, 5)], dtype=np.double)
  108. points_rescaled = np.array([(0,0), (0,1), (1,1), (1,0), (0.1, 0.05)], dtype=np.double)
  109. values = np.array([1., 2., -3., 5., 9.], dtype=np.double)
  110. xx, yy = np.broadcast_arrays(np.linspace(0, 10, 14)[:,None],
  111. np.linspace(0, 100, 14)[None,:])
  112. xx = xx.ravel()
  113. yy = yy.ravel()
  114. xi = np.array([xx, yy]).T.copy()
  115. for method in ('nearest', 'linear', 'cubic'):
  116. msg = method
  117. zi = griddata(points_rescaled, values, xi/np.array([10, 100.]),
  118. method=method)
  119. zi_rescaled = griddata(points, values, xi, method=method,
  120. rescale=True)
  121. assert_allclose(zi, zi_rescaled, err_msg=msg,
  122. atol=1e-12)
  123. def test_xi_1d(self):
  124. # Check that 1-D xi is interpreted as a coordinate
  125. x = np.array([(0,0), (-0.5,-0.5), (-0.5,0.5), (0.5, 0.5), (0.25, 0.3)],
  126. dtype=np.double)
  127. y = np.arange(x.shape[0], dtype=np.double)
  128. y = y - 2j*y[::-1]
  129. xi = np.array([0.5, 0.5])
  130. for method in ('nearest', 'linear', 'cubic'):
  131. p1 = griddata(x, y, xi, method=method)
  132. p2 = griddata(x, y, xi[None,:], method=method)
  133. assert_allclose(p1, p2, err_msg=method)
  134. xi1 = np.array([0.5])
  135. xi3 = np.array([0.5, 0.5, 0.5])
  136. assert_raises(ValueError, griddata, x, y, xi1,
  137. method=method)
  138. assert_raises(ValueError, griddata, x, y, xi3,
  139. method=method)
  140. class TestNearestNDInterpolator:
  141. def test_nearest_options(self):
  142. # smoke test that NearestNDInterpolator accept cKDTree options
  143. npts, nd = 4, 3
  144. x = np.arange(npts*nd).reshape((npts, nd))
  145. y = np.arange(npts)
  146. nndi = NearestNDInterpolator(x, y)
  147. opts = {'balanced_tree': False, 'compact_nodes': False}
  148. nndi_o = NearestNDInterpolator(x, y, tree_options=opts)
  149. assert_allclose(nndi(x), nndi_o(x), atol=1e-14)
  150. def test_nearest_list_argument(self):
  151. nd = np.array([[0, 0, 0, 0, 1, 0, 1],
  152. [0, 0, 0, 0, 0, 1, 1],
  153. [0, 0, 0, 0, 1, 1, 2]])
  154. d = nd[:, 3:]
  155. # z is np.array
  156. NI = NearestNDInterpolator((d[0], d[1]), d[2])
  157. assert_array_equal(NI([0.1, 0.9], [0.1, 0.9]), [0, 2])
  158. # z is list
  159. NI = NearestNDInterpolator((d[0], d[1]), list(d[2]))
  160. assert_array_equal(NI([0.1, 0.9], [0.1, 0.9]), [0, 2])
  161. class TestNDInterpolators:
  162. @parametrize_interpolators
  163. def test_broadcastable_input(self, interpolator):
  164. # input data
  165. np.random.seed(0)
  166. x = np.random.random(10)
  167. y = np.random.random(10)
  168. z = np.hypot(x, y)
  169. # x-y grid for interpolation
  170. X = np.linspace(min(x), max(x))
  171. Y = np.linspace(min(y), max(y))
  172. X, Y = np.meshgrid(X, Y)
  173. XY = np.vstack((X.ravel(), Y.ravel())).T
  174. interp = interpolator(list(zip(x, y)), z)
  175. # single array input
  176. interp_points0 = interp(XY)
  177. # tuple input
  178. interp_points1 = interp((X, Y))
  179. interp_points2 = interp((X, 0.0))
  180. # broadcastable input
  181. interp_points3 = interp(X, Y)
  182. interp_points4 = interp(X, 0.0)
  183. assert_equal(interp_points0.size ==
  184. interp_points1.size ==
  185. interp_points2.size ==
  186. interp_points3.size ==
  187. interp_points4.size, True)
  188. @parametrize_interpolators
  189. def test_read_only(self, interpolator):
  190. # input data
  191. np.random.seed(0)
  192. xy = np.random.random((10, 2))
  193. x, y = xy[:, 0], xy[:, 1]
  194. z = np.hypot(x, y)
  195. # interpolation points
  196. XY = np.random.random((50, 2))
  197. xy.setflags(write=False)
  198. z.setflags(write=False)
  199. XY.setflags(write=False)
  200. interp = interpolator(xy, z)
  201. interp(XY)