Implement limit to bounds in GridInterpolator

py_wake/deficit_models/fuga.py

@@ -81,8 +81,9 @@ class FugaUtils():
         return np.concatenate([((1, -1)[anti_symmetric]) * x[::-1], x[1:]])
 
     def load_luts(self, UVLT=['UL', 'UT', 'VL', 'VT'], zlevels=None):
-        return np.array([[np.fromfile(str(self.path / (self.prefix + '%04d%s.dat' % (j, uvlt))), np.dtype('<f'), -1)
-                          for j in (zlevels or self.zlevels)] for uvlt in UVLT])
+        luts = np.array([[np.fromfile(str(self.path / (self.prefix + '%04d%s.dat' % (j, uvlt))), np.dtype('<f'), -1)
+                          for j in (zlevels or self.zlevels)] for uvlt in UVLT]).astype(np.float)
+        return luts.reshape((len(UVLT), len(zlevels or self.zlevels), self.ny // 2, self.nx))
 
 
 class FugaDeficit(DeficitModel, FugaUtils):

py_wake/site/xrsite.py

@@ -10,12 +10,14 @@ class XRSite(Site):
     use_WS_bins = False
 
     def __init__(self, ds, initial_position=None, interp_method='linear', shear=None, distance=StraightDistance(),
-                 default_ws=np.arange(3, 26)):
-        if interp_method not in ['linear', 'nearest']:
-            raise NotImplementedError('interp_method=%s not implemented' % interp_method)
+                 default_ws=np.arange(3, 26), bounds='check'):
+        assert interp_method in [
+            'linear', 'nearest'], 'interp_method "%s" not implemented. Must be "linear" or "nearest"' % interp_method
+        assert bounds in ['check', 'limit', 'ignore'], 'bounds must be "check", "limit" or "ignore"'
 
         self.interp_method = interp_method
         self.shear = shear
+        self.bounds = bounds
 
         Site.__init__(self, distance)
         self.default_ws = default_ws
@@ -127,7 +129,7 @@ class XRSite(Site):
 
         if len(ip_dims) > 0:
             grid_interp = GridInterpolator([var.coords[k].data for k in ip_dims], data,
-                                           method=self.interp_method)
+                                           method=self.interp_method, bounds=self.bounds)
 
             # get dimension of interpolation coordinates
             I = (1, len(coords.get('x', coords.get('y', coords.get('h', coords.get('i', [None]))))))[

py_wake/tests/test_sites/test_site.py

@@ -57,7 +57,7 @@ def test_elevation(site):
 
 
 def test_missing_interp_method():
-    with pytest.raises(NotImplementedError, match="interp_method=missing_method not implemented"):
+    with pytest.raises(AssertionError, match='interp_method "missing_method" not implemented. Must be "linear" or "nearest"'):
         site = UniformWeibullSite([1], [10], [2], .75, interp_method='missing_method')
 
 

py_wake/tests/test_utils/test_grid_interpolator.py

@@ -15,11 +15,11 @@ def test_grid_interpolator_not_match():
 
 
 def test_grid_interpolator_wrong_method():
-    with pytest.raises(ValueError, match='Method must be "linear" or "nearest"'):
+    with pytest.raises(AssertionError, match='method must be "linear" or "nearest"'):
         GridInterpolator([[10, 20, 30]], [1, 2, 3], method="cubic")(None)
 
 
-def test_grid_interpolator_outside_area():
+def test_grid_interpolator_outside_area_bounds_check():
     with pytest.raises(ValueError, match='Point 0, dimension 0 with value 9.900000 is outside range 10.000000-30.000000'):
         GridInterpolator([[10, 20, 30]], [1, 2, 3])(9.9)
 
@@ -36,6 +36,16 @@ def test_grid_interpolator_outside_area():
         GridInterpolator([[5, 10, 15], [200, 300]], np.arange(6).reshape(3, 2))([(5, 199.9), (4.8, 200), (5, 301)])
 
 
+def test_grid_interpolator_outside_area_bounds_limit():
+    gi = GridInterpolator([[10, 20, 30]], [1, 2, 3])
+    npt.assert_array_equal(gi([[10]]), gi([[9.9]], bounds='limit'))
+    npt.assert_array_equal(gi([[30]]), gi([[30.1]], bounds='limit'))
+    npt.assert_array_equal(gi([[20], [30]]), gi([[20], [30.1]], bounds='limit'))
+    gi = GridInterpolator([[5, 10, 15], [200, 300]], np.arange(6).reshape(3, 2))
+    npt.assert_array_equal(gi([(5, 300)]), gi([(5, 300.1)], bounds='limit'))
+    npt.assert_array_equal(gi([(5, 200), (5, 200), (5, 300)]), gi([(5, 199.9), (4.8, 200), (5, 301)], bounds='limit'))
+
+
 def test_grid_interpolator_2d():
     x = [5, 10, 15]
     y = [200, 300, 400, 500]

py_wake/utils/grid_interpolator.py

@@ -5,8 +5,26 @@ from py_wake.tests import npt
 
 class GridInterpolator(object):
     # Faster than scipy.interpolate.interpolate.RegularGridInterpolator
-    def __init__(self, x, V, method='linear'):
+    def __init__(self, x, V, method='linear', bounds='check'):
+        """
+        Parameters
+        ----------
+        x : list of array_like
+            Interpolation coordinates
+        V : array_like
+            Interpolation data (more dimensions than coordinates is possible)
+        method : {'linear' or 'nearest} or None
+            Overrides self.method
+        bounds : {'check', 'limit', 'ignore'}
+            Specifies how bounds is handled:\n
+            - 'check': bounds check is performed. An error is raised if interpolation point outside area
+            - 'limit': interpolation points are forced inside the area
+            - 'ignore': Faster option with no check. Use this option if data is guaranteed to be inside the area
+        """
         self.x = x
+        self.V = V
+        self.bounds = bounds
+        self.method = method
         self.n = np.array([len(x) for x in x])
         self.x0 = np.array([x[0] for x in x])
         dx = np.array([x[min(1, len(x) - 1)] - x[0] for x in x])
@@ -24,12 +42,22 @@ class GridInterpolator(object):
         ui[:, dx == 0] = 0
         self.ui = ui
 
-        self.method = method
+    def __call__(self, xp, method=None, bounds=None):
+        """Interpolate points
 
-    def __call__(self, xp, method=None, check_bounds=True):
+        Parameters
+        ----------
+        xp : array_like
+            Interpolation points, shape=(n_points, interpolation_dimensions)
+        method : {'linear' or 'nearest} or None
+            Overrides self.method if not None
+        bounds : {'check', 'limit', 'ignore'} or None
+            Overrides self.bounds if not None
+        """
         method = method or self.method
-        if method not in ['linear', 'nearest']:
-            raise ValueError('Method must be "linear" or "nearest"')
+        bounds = bounds or self.bounds
+        assert method in ['linear', 'nearest'], 'method must be "linear" or "nearest"'
+        assert bounds in ['check', 'limit', 'ignore'], 'bounds must be "check", "limit" or "ignore"'
         xp = np.asarray(xp)
         xpi = (xp - self.x0) / self.dx
         if len(self.irregular_axes):
@@ -41,13 +69,15 @@ class GridInterpolator(object):
             irreg_dx = irreg_x1 - irreg_x0
             xpi[:, self.irregular_axes] = irreg_i.T + (xp[:, self.irregular_axes] - irreg_x0.T) / irreg_dx.T
 
-        if check_bounds and (np.any(xpi < 0) or np.any(xpi + 1 > self.n[na])):
+        if bounds == 'check' and (np.any(xpi < 0) or np.any(xpi + 1 > self.n[na])):
             if -xpi.min() > (xpi + 1 - self.n[na]).max():
                 point, dimension = np.unravel_index(xpi.argmin(), np.atleast_2d(xpi).shape)
             else:
                 point, dimension = np.unravel_index(((xpi + 1 - self.n[na])).argmax(), np.atleast_2d(xpi).shape)
             raise ValueError("Point %d, dimension %d with value %f is outside range %f-%f" %
                              (point, dimension, np.atleast_2d(xp)[point, dimension], self.x[dimension][0], self.x[dimension][-1]))
+        if bounds == 'limit':
+            xpi = np.minimum(np.maximum(xpi, 0), self.n - 1)
         xpi0 = xpi.astype(int)
         xpif = xpi - xpi0
         if method == 'nearest':