Add Nearest-Neighbor interpolation to heterogeneous wind map

floris/core/flow_field.py

@@ -2,10 +2,11 @@
 from __future__ import annotations
 
 import attrs
+import copy
 import matplotlib.path as mpltPath
 import numpy as np
 from attrs import define, field
-from scipy.interpolate import LinearNDInterpolator
+from scipy.interpolate import LinearNDInterpolator, NearestNDInterpolator
 from scipy.spatial import ConvexHull
 from shapely.geometry import Polygon
 
@@ -286,30 +287,45 @@ def generate_heterogeneous_wind_map(self):
         y = self.heterogeneous_inflow_config['y']
         z = self.heterogeneous_inflow_config['z']
 
+        if "interp_method" in self.heterogeneous_inflow_config.keys():
+            interp_method = self.heterogeneous_inflow_config["interp_method"]
+        else:
+            interp_method = "linear"
+
         if z is not None:
             # Compute the 3-dimensional interpolants for each wind direction
             # Linear interpolation is used for points within the user-defined area of values,
             # while the freestream wind speed is used for points outside that region
-            in_region = [
-                self.interpolate_multiplier_xyz(x, y, z, multiplier, fill_value=1.0)
-                for multiplier in speed_multipliers
-            ]
+            F = self.interpolate_multiplier_xyz(x, y, z, speed_multipliers[0], fill_value=1.0, interp_method=interp_method)
+            in_region = []
+            for multiplier in speed_multipliers:
+                if interp_method == "linear":
+                    F.values = multiplier[:, None]  # Linear
+                else: 
+                    F.values = multiplier  # NN
+                in_region.append(copy.deepcopy(F))
+
         else:
             # Compute the 2-dimensional interpolants for each wind direction
             # Linear interpolation is used for points within the user-defined area of values,
             # while the freestream wind speed is used for points outside that region
-            in_region = [
-                self.interpolate_multiplier_xy(x, y, multiplier, fill_value=1.0)
-                for multiplier in speed_multipliers
-            ]
+            self.interpolate_multiplier_xy(x, y, multiplier, fill_value=1.0, interp_method=interp_method)
+            in_region = []
+            for multiplier in speed_multipliers:
+                if interp_method == "linear":
+                    F.values = multiplier[:, None]  # Linear
+                else: 
+                    F.values = multiplier  # NN
+                in_region.append(copy.deepcopy(F))
 
         self.het_map = in_region
 
     @staticmethod
     def interpolate_multiplier_xy(x: NDArrayFloat,
                                   y: NDArrayFloat,
                                   multiplier: NDArrayFloat,
-                                  fill_value: float = 1.0):
+                                  fill_value: float = 1.0,
+                                  interp_method: str = "linear"):
         """Return an interpolant for a 2D multiplier field.
 
         Args:
@@ -321,16 +337,21 @@ def interpolate_multiplier_xy(x: NDArrayFloat,
         Returns:
             LinearNDInterpolator: interpolant
         """
-
-        return LinearNDInterpolator(list(zip(x, y)), multiplier, fill_value=fill_value)
+        if interp_method == "linear":
+            return LinearNDInterpolator(list(zip(x, y)), multiplier, fill_value=fill_value)
+        elif interp_method == "nn":
+            return NearestNDInterpolator(list(zip(x, y)), multiplier)
+        else:
+            raise UserWarning("Incompatible interpolation method specified.")
 
 
     @staticmethod
     def interpolate_multiplier_xyz(x: NDArrayFloat,
                                    y: NDArrayFloat,
                                    z: NDArrayFloat,
                                    multiplier: NDArrayFloat,
-                                   fill_value: float = 1.0):
+                                   fill_value: float = 1.0,
+                                   interp_method: str = "linear"):
         """Return an interpolant for a 3D multiplier field.
 
         Args:
@@ -344,4 +365,9 @@ def interpolate_multiplier_xyz(x: NDArrayFloat,
             LinearNDInterpolator: interpolant
         """
 
-        return LinearNDInterpolator(list(zip(x, y, z)), multiplier, fill_value=fill_value)
+        if interp_method == "linear":
+            return LinearNDInterpolator(list(zip(x, y, z)), multiplier, fill_value=fill_value)
+        elif interp_method == "nn":
+            return NearestNDInterpolator(list(zip(x, y, z)), multiplier)
+        else:
+            raise UserWarning("Incompatible interpolation method specified.")

floris/heterogeneous_map.py

@@ -29,7 +29,11 @@ class HeterogeneousMap(LoggingManager):
             (degrees). Optional.
         wind_speeds (NDArrayFloat, optional): A 1D NumPy array (size num_ws) of wind speeds (m/s).
             Optional.
-
+    	interp_method (str, optional): Interpolation method used to calculate the heterogeneous
+            wind speeds at various locations in the wind farm. Options are 'linear' and 'nn',
+            representing linear interpolation and nearest-neighbor interpolation respectively.
+            Linear interpolation is accurate, nearest-neighbor interpolation is very fast but
+            inaccurate. Defaults to 'linear'. Optional.
 
     Notes:
         * If wind_directions and wind_speeds are both defined, then they must be the same length
@@ -45,6 +49,7 @@ def __init__(
         z: NDArrayFloat = None,
         wind_directions: NDArrayFloat = None,
         wind_speeds: NDArrayFloat = None,
+        interp_method: str = "linear",
     ):
         # Check that x, y and speed_multipliers are lists or numpy arrays
         if not isinstance(x, (list, np.ndarray)):
@@ -62,6 +67,7 @@ def __init__(
         self.x = np.array(x)
         self.y = np.array(y)
         self.speed_multipliers = np.array(speed_multipliers)
+        self.interp_method = str(interp_method)
 
         # If z is provided, save it as an np array
         if z is not None:
@@ -251,13 +257,15 @@ def get_heterogeneous_inflow_config(
                 "x": self.x,
                 "y": self.y,
                 "speed_multipliers": speed_multipliers_by_findex,
+                "interp_method": self.interp_method,
             }
         else:
             return {
                 "x": self.x,
                 "y": self.y,
                 "z": self.z,
                 "speed_multipliers": speed_multipliers_by_findex,
+                "interp_method": self.interp_method,
             }
 
     def get_heterogeneous_map_2d(self, z: float):
@@ -287,6 +295,7 @@ def get_heterogeneous_map_2d(self, z: float):
             speed_multipliers=speed_multipliers,
             wind_directions=self.wind_directions,
             wind_speeds=self.wind_speeds,
+            interp_method=self.interp_method,
         )
 
     @staticmethod