ISO noise model provided by Camilla Nyborg implemented and integrated into simulation result

py_wake/examples/data/swt_dd_142_4100_noise/swt_dd_142_4100.py

+import pandas as pd
+from py_wake.wind_turbines.power_ct_functions import PowerCtXr
+from py_wake.wind_turbines._wind_turbines import WindTurbine
+import xarray as xr
+import numpy as np
+from py_wake.examples.data import swt_dd_142_4100_noise
+import os
+from py_wake.utils.grid_interpolator import GridInterpolator
+from py_wake.utils.plotting import setup_plot
+
+
+class SWT_DD_142_4100(WindTurbine):
+    """Siemens SWT-DD-142, 4.1MW, including sound power levels
+
+    Data extracted from Windpro and stored as netcdf in SWT-DD-142_4100.nc
+    """
+
+    def __init__(self):
+        ds = xr.load_dataset(os.path.dirname(swt_dd_142_4100_noise.__file__) + '/SWT-DD-142_4100.nc')
+        power_ct = PowerCtXr(ds, 'kW')
+        ip = GridInterpolator([ds.mode.values, ds.ws.values], ds.SoundPower.transpose('mode', 'ws', 'freq').values,
+                              method=['nearest', 'linear'])
+
+        def sound_power_level(ws, mode, **_):
+            mode = np.zeros_like(ws) + mode
+            return ds.freq.values, ip(np.array([np.round(mode).flatten(), ws.flatten()]).T).reshape(ws.shape + (-1,))
+
+        WindTurbine.__init__(self, name='SWT-DD-142', diameter=142, hub_height=109,
+                             powerCtFunction=power_ct, sound_power_level=sound_power_level)
+
+
+def make_netcdf():
+    # temp function used to convert the data from excel to netcdf
+
+    # Read in the turbine power, ct and noise curves retrieved from Windpro
+    SWT = pd.read_excel(r'../../../../SWT-DD-142_4100.xlsx')  # Siemens turbine with different modes
+    modes = np.arange(7)
+    ws = SWT[SWT.Mode == 0]['WindSpeed'].values.astype(float)
+    freqs = [63, 125, 250, 500, 1000, 2000, 4000, 8000]
+
+    SWT = SWT.rename(columns={'f8000Hz ': 'f8000Hz'})
+
+    ds = xr.Dataset({'SoundPower': (('mode', 'freq', 'ws'),
+                                    [[SWT[SWT.Mode == mode][f'f{freq}Hz'].values.astype(float) for freq in freqs]
+                                     for mode in modes]),
+                     **{k: (('mode', 'ws'),
+                            [SWT[SWT.Mode == mode][k].values.astype(float) for mode in modes])
+                        for k in ['LwaRef', 'Power', 'Ct']}},
+                    coords={'ws': ws, 'freq': freqs, 'mode': modes})
+    ds.to_netcdf('SWT-DD-142_4100.nc')
+
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    wt = SWT_DD_142_4100()
+
+    ws = np.arange(3, 26)
+    for m in range(7):
+        plt.plot(ws, wt.power(ws, mode=m) / 1000, label=f'mode: {m}')
+    setup_plot(xlabel='Wind speed [m/s]', ylabel='Power [kW]')
+
+    plt.figure()
+    for m in range(7):
+        freq, sound_power = wt.sound_power(ws=10, mode=m)
+        plt.plot(freq, sound_power[0], label=f'mode: {m}')
+    setup_plot(xlabel='Frequency [Hz]', ylabel='Sound power [dB]')
+    plt.show()

py_wake/noise_models/iso.py

+import numpy as np
+from py_wake.site.distance import StraightDistance
+from numpy import newaxis as na
+from py_wake.utils.grid_interpolator import GridInterpolator
+
+
+class ISONoiseModel:
+
+    def __init__(self, src_x, src_y, src_h, freqs, sound_power_level, elevation_function=None):
+        """ISONoise model based on
+
+        DSF/ISO/DIS 9613-2
+        Acoustics – Attenuation of sound during propagation – Part 2:
+        Engineering method for the prediction of sound pressure levels outdoors
+
+        DS/ISO 9613-1:1993
+        Akustik. Måling og beskrivelse af ekstern støj. Lydudbredelsesdæmpning udendørs. Del 1:
+        Metode til beregning af luftabsorption
+
+        The code is a vectorized version of the implementaion made by Camilla Marie Nyborg <cmny@dtu.dk>
+
+        The model models the sound pressure level from a number of sound sources at a number of receivers taking into
+        account
+        - spherical geometrical spreading
+        - ground reflection/absorption
+        - atmospheric absorption (DS/ISO 9613-1:1993)
+
+
+        Parameters
+        ----------
+        src_x : array_like
+            x coordinate of sound sources
+        src_y : array_like
+            y coordinate of sound sources
+        src_h : array_like or float
+            height of sound sources (typically wind turbine hub height)
+        freqs : array_like
+            Frequencies of sound_power_level
+        sound_power_level : array_like
+            Emitted sound power level, dim=(n_sources, n_freq)
+
+        """
+
+        if elevation_function:
+            src_h += elevation_function(src_x, src_y)
+        self.elevation_function = elevation_function
+        self.src_x, self.src_y = src_x, src_y
+        self.src_h = np.zeros_like(src_x) + src_h
+        self.n_src = len(src_x)
+        self.distance = StraightDistance()
+        self.freqs = freqs
+        self.sound_power_level = sound_power_level
+
+    def ground_eff(self, ground_distance_ij, distance_ij, ground_type):
+        # Ground effects ISO
+
+        freqs_init = np.array([63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0, 8000.0])
+
+        src_h_ij = self.src_h[:, na]
+        rec_h_ij = self.distance.dst_h_j[na]
+
+        hei_check_ij = 30.0 * (src_h_ij + rec_h_ij)
+
+        q_ij = np.where(ground_distance_ij <= hei_check_ij, 0, 1 - hei_check_ij / ground_distance_ij)
+
+        Gs = Gr = Gm = ground_type
+
+        # for f=63Hz Am =-3*q, for other frequencies Am = -3 * q * (1-Gm)
+        fak = np.where(freqs_init == 63.0, 1, (1 - Gm))
+        Am_fij = -3.0 * q_ij[na] * fak[:, na, na]
+
+        # 125 Hz
+        a_source_ij = 1.5 + 3.0 * np.exp(-0.12 * (src_h_ij - 5.0)**2) * (1.0 - np.exp(-ground_distance_ij / 50.0)) + \
+            5.7 * np.exp(-0.09 * src_h_ij**2) * \
+            (1.0 - np.exp(-2.8 * (10.0**(-6.0)) * (ground_distance_ij**2)))
+        a_rec_ij = 1.5 + 3.0 * np.exp(-0.12 * (rec_h_ij - 5.0)**2) * (1.0 - np.exp(-ground_distance_ij / 50.0)) + \
+            5.7 * np.exp(-0.09 * rec_h_ij**2) * (1.0 - np.exp(-2.8 * (10.0**(-6.0)) * (ground_distance_ij**2)))
+
+        # 250 Hz
+        b_source_ij = 1.5 + 8.6 * np.exp(-0.09 * src_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
+        b_rec_ij = 1.5 + 8.6 * np.exp(-0.09 * rec_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
+
+        # 500 Hz
+        c_source_ij = 1.5 + 14.0 * np.exp(-0.46 * src_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
+        c_rec_ij = 1.5 + 14.0 * np.exp(-0.46 * rec_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
+
+        # 1000 Hz
+        d_source_ij = 1.5 + 5.0 * np.exp(-0.9 * src_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
+        d_rec_ij = 1.5 + 5.0 * np.exp(-0.9 * rec_h_ij**2) * (1.0 - np.exp(-ground_distance_ij / 50.0))
+
+        zeros = np.zeros_like(ground_distance_ij)
+        As_fij = np.array([zeros - 1.5,  # 63Hz
+                           -1.5 + Gs * a_source_ij,  # 125Hz
+                           -1.5 + Gs * b_source_ij,  # 250Hz
+                           -1.5 + Gs * c_source_ij,  # 500Hz
+                           -1.5 + Gs * d_source_ij,  # 1000Hz
+                           zeros - 1.5 * (1.0 - Gs),  # 2000Hz
+                           zeros - 1.5 * (1.0 - Gs)  # 8000Hz
+                           ])
+
+        Ar_fij = np.array([zeros - 1.5,  # 63 Hz
+                           -1.5 + Gr * a_rec_ij,  # 125 Hz
+                           -1.5 + Gr * b_rec_ij,  # 250 Hz
+                           -1.5 + Gr * c_rec_ij,  # 500 Hz
+                           -1.5 + Gr * d_rec_ij,  # 1000 Hz
+                           zeros - 1.5 * (1.0 - Gr),  # 2000 Hz
+                           zeros - 1.5 * (1.0 - Gr)  # 8000 Hz
+                           ])
+
+        # Interpolation to other frequencies are actually not following the standard completely.
+        ip = GridInterpolator([freqs_init], np.moveaxis([As_fij, Ar_fij, Am_fij], 0, -1))
+
+        # interpolate to freq and sum up As+Ar+Am
+        Agr_ijf = np.moveaxis(ip(np.array([self.freqs]).T).sum(-1), 0, -1)
+
+        # Area of sphere: 4pi R^2
+        # 10*log_10(4pi) ~ 11
+        # 10 * log_10(R^2) = 20 * log_10(R) =
+        # reference area = 1.0 m^2
+        Adiv_ij = 20.0 * np.log10(distance_ij / 1.0) + 11.0  # The geometrical spreading (divergence)
+
+        ISO_ground_ijf = - Adiv_ij[:, :, na] - Agr_ijf
+
+        return ISO_ground_ijf
+
+    def atmab(self, distance_ij, T0, RH0):
+        # Atmospheric absorption
+        T_0 = 293.15
+        T_01 = 273.16
+        T = T0 + 273.15
+        p_s0 = 1.0
+
+        psat = p_s0 * 10.0**(-6.8346 * (T_01 / T)**1.261 + 4.6151)
+        h = p_s0 * (RH0) * (psat / p_s0)
+        F_rO = 1.0 / p_s0 * (24.0 + 4.04 * (10**4.0) * h * (0.02 + h) / (0.391 + h))
+        F_rN = 1.0 / p_s0 * (T_0 / T)**(0.5) * (9.0 + 280 * h *
+                                                np.exp(-4.17 * ((T_0 / T)**(1.0 / 3.0) - 1.0)))
+        alpha_ps = self.freqs**2.0 / p_s0 * (1.84 * (10**(-11)) * (T / T_0)**(0.5) + (T / T_0)**(-5.0 / 2.0) *
+                                             (0.01275 * np.exp(-2239.1 / T) / (F_rO + self.freqs**2.0 / F_rO) +
+                                              0.1068 * np.exp(-3352 / T) / (F_rN + self.freqs**2.0 / F_rN)))
+
+        ISO_alpha = alpha_ps[na, na] * 20.0 / np.log(10.0) * distance_ij[:, :, na]
+        return ISO_alpha
+
+    def transmission_loss(self, rec_x, rec_y, rec_h, ground_type, Temp, RHum):
+        # transmission loss = ground effects + atmospheric absorption
+        rec_h = np.zeros_like(rec_x) + rec_h
+        if self.elevation_function:
+            rec_h += self.elevation_function(rec_x, rec_y)
+        self.distance.setup(self.src_x, self.src_y, self.src_h, [rec_x, rec_y, rec_h])
+        ground_distance_ij = np.sqrt(self.distance.dx_ij**2 + self.distance.dy_ij**2)
+        distance_ij = np.sqrt(self.distance.dx_ij**2 + self.distance.dy_ij**2 + self.distance.dh_ij**2)
+
+        atm_abs_ijf = self.atmab(distance_ij, T0=Temp, RH0=RHum)  # The atmospheric absorption term
+        ground_eff_ijf = self.ground_eff(ground_distance_ij, distance_ij, ground_type)
+
+        return ground_eff_ijf - atm_abs_ijf  # Delta_SPL
+
+    def __call__(self, rec_x, rec_y, rec_h, Temp, RHum, ground_type):
+        """Calculate the sound pressure level at a list of reveicers
+
+        Parameters
+        ----------
+        rec_x : array_like
+            x coordinate, [m], of receivers
+        rec_y : array_like
+            y coordinate, [m], of receivers
+        rec_h : array_like or float
+            height, [m], of receivers (typically 2m)
+        Temp : float
+            Temperature (deg celcius)
+        RHum : float
+            Relative humidity (%)
+        ground_type : float
+            Factor describing the amount of ground absorption, 0=hard reflective ground,
+            e.g. paving, water, ice and concrete. 1= soft porous absorbing ground,
+            e.g. grass, trees, other vegetation and farming land
+
+        Returns
+        -------
+        total_spl : array_like
+            Total sound pressure level at each receiver
+        spl : array_like
+            Sound pressure level of freqs, dim=(n_receiver, n_freq)
+        """
+
+        # Computing total transmission loss
+        Delta_SPL_ijf = self.transmission_loss(rec_x, rec_y, rec_h, ground_type, Temp, RHum)
+        sound_power_ijxxf = self.sound_power_level[:, na]
+        I, J, F = Delta_SPL_ijf.shape
+        shape = [I, J] + [1] * (len(sound_power_ijxxf.shape) - len(Delta_SPL_ijf.shape)) + [F]
+        Delta_SPL_ijxxf = Delta_SPL_ijf.reshape(shape)
+
+        # Add negative transmission loss to emitted sound and sum over sources to get sound pressure level
+        spl_jxxf = 10.0 * np.log10(np.sum(10.0**((self.sound_power_level[:, na] + Delta_SPL_ijxxf) / 10.0), axis=0))
+
+        # Sum over frequencies to get total sound pressure level
+        total_spl_jxx = 10.0 * np.log10(np.sum(10.0**(spl_jxxf / 10.0), axis=-1))
+        return total_spl_jxx, spl_jxxf
+
+
+def main():
+    if __name__ == '__main__':
+        import matplotlib.pyplot as plt
+        from py_wake.deficit_models.gaussian import ZongGaussian
+        from py_wake.flow_map import XYGrid
+        from py_wake.turbulence_models.crespo import CrespoHernandez
+
+        from py_wake.site._site import UniformSite
+        from py_wake.examples.data.swt_dd_142_4100_noise.swt_dd_142_4100 import SWT_DD_142_4100
+        from py_wake.utils.layouts import rectangle
+        from py_wake.utils.plotting import setup_plot
+
+        wt = SWT_DD_142_4100()
+        wfm = ZongGaussian(UniformSite(), wt, turbulenceModel=CrespoHernandez())
+        x, y = rectangle(5, 5, 5 * wt.diameter())
+        sim_res = wfm(x, y, wd=270, ws=8, mode=0)
+        nm = sim_res.noise_model()
+        ax1, ax2, ax3 = plt.subplots(1, 3, figsize=(16, 4))[1]
+        sim_res.flow_map().plot_wake_map(ax=ax1)
+
+        ax1.plot([x[0]], [1000], '.', label='Receiver 1')
+        ax1.plot([x[-1]], [1000], '.', label='Receiver 2')
+        ax1.legend()
+        total_sp_jlk, spl_jlkf = nm(rec_x=[x[0], x[-1]], rec_y=[1000, 1000], rec_h=2, Temp=20, RHum=80, ground_type=0.0)
+        ax2.plot(nm.freqs, spl_jlkf[0, 0, 0], label='Receiver 1')
+        ax2.plot(nm.freqs, spl_jlkf[1, 0, 0], label='Receiver 2')
+        setup_plot(xlabel='Frequency [Hz]', ylabel='Sound pressure level [dB]', ax=ax2)
+        plt.tight_layout()
+
+        nmap = sim_res.noise_map(grid=XYGrid(x=np.linspace(-1000, 5000, 100), y=np.linspace(-1000, 1000, 50), h=2))
+        nmap['Total sound pressure level'].squeeze().plot(ax=ax3)
+        wt.plot(x, y, ax=ax3)
+
+        plt.show()
+
+
+main()

py_wake/site/xrsite.py

@@ -124,7 +124,7 @@ class XRSite(Site):
                                                np.asarray(y_i, dtype=(float, np.complex128)[np.iscomplexobj(y_i)]),
                                                mode='valid')
         else:
-            return x_i * 0
+            return np.asarray(x_i) * 0
 
     def _local_wind(self, localWind, ws_bins=None):
         """

py_wake/tests/test_noise_models/test_iso.py

+import numpy as np
+from py_wake.noise_models.iso import ISONoiseModel
+from py_wake.tests import npt
+from numpy import newaxis as na
+from py_wake.deficit_models.gaussian import BastankhahGaussian
+from py_wake.examples.data.swt_dd_142_4100_noise.swt_dd_142_4100 import SWT_DD_142_4100
+from py_wake.site._site import UniformSite
+from py_wake.utils.layouts import rectangle
+from py_wake.utils.plotting import setup_plot
+from py_wake.flow_map import XYGrid
+
+
+def test_iso_noise_model():
+
+    Temp = 20  # Temperature
+    RHum = 80.0  # Relative humidity
+
+    rec_x = [2000, 3000, 4000]  # receiver xy-position
+    rec_y = [0, 0, 0]
+
+    source_x = [0, 0]  # source xy-position
+    source_y = [0, 500]
+
+    freqs = np.array([63.0, 125.0, 250.0, 500.0, 1000.0, 2000.0, 4000.0, 8000.0])
+
+    # If the source and receiver are placed in complex
+    # terrain then the height of the terrain should be added to these heights.
+    z_hub = 100  # source height
+    z_rec = 2  # receiver height
+
+    ground_type = 1  # Ranging from 0 to 1 (hard to soft ground)
+
+    sound_power_level = np.array([[89.4, 93.6, 97.2, 98.6, 101., 102.3, 96.7, 84.1],
+                                  [89.2, 93.2, 96.2, 97.6, 100., 101.3, 95.7, 83.1]])
+    iso = ISONoiseModel(src_x=source_x, src_y=source_y, src_h=z_hub, freqs=freqs, sound_power_level=sound_power_level)
+
+    Delta_SPL = iso.transmission_loss(rec_x, rec_y, z_rec, ground_type, Temp, RHum)
+    delta_SPL_ref = [[[-74.18868997998351, -82.6237111823142, -85.106899362965, -84.78075941279131, -87.47936788035022,
+                       -95.0531580482448, -119.90461076531315, -216.18076560019855],
+                      [-77.78341234414849, -86.43781608834009, -89.6588031834609, -91.0544347291193, -96.14098255652445,
+                       -107.56227940740277, -144.8146479691885, -289.1327626666617],
+                      [-79.65387282494626, -89.23269240205617, -93.19182772534491, -96.30991899366798, -103.78536068849834,
+                       -119.05570214846978, -168.71394345143312, -361.09322135290586]],
+                     [[-74.4562086405804, -82.90475257479693, -85.43331381258002, -85.21311520726341, -88.05865461649006,
+                       -95.86918353475227, -121.48366996751929, -220.71586737235828],
+                      [-77.90553623031623, -86.56901853514375, -89.82054728871552, -91.28744741509202, -96.47283823614013,
+                         -108.0533933516973, -145.81906793231144, -292.12576375254764],
+                      [-79.70589475152278, -89.3092674680563, -93.29138286495427, -96.46309784962982, -104.01291072740128,
+                         -119.40308086820056, -169.44754252350324, -363.32306335554176]]]
+    npt.assert_array_almost_equal(delta_SPL_ref, Delta_SPL, 8)
+    total_spl, spl = iso(rec_x, rec_y, z_rec, ground_type=ground_type, Temp=Temp, RHum=RHum)
+    npt.assert_array_almost_equal([23.068816073636583, 18.034141554900494, 15.043308370722233],
+                                  total_spl)
+
+
+def test_iso_noise_map():
+    wt = SWT_DD_142_4100()
+    wfm = BastankhahGaussian(UniformSite(), wt)
+    x, y = rectangle(5, 5, 5 * wt.diameter())
+    sim_res = wfm(x, y, wd=270, ws=8, mode=0)
+
+    nm = sim_res.noise_model()
+    total_spl_jlk, spl_jlkf = nm(rec_x=[x[0], x[-1]], rec_y=[1000, 1000], rec_h=2, Temp=20, RHum=80, ground_type=0.0)
+    sim_res.noise_map()  # cover default grid
+    nmap = sim_res.noise_map(grid=XYGrid(x=np.linspace(-1000, 4000, 100), y=np.linspace(-1000, 1000, 50), h=2))
+
+    npt.assert_array_almost_equal(total_spl_jlk.squeeze(), [34.87997084, 29.56191044])
+    npt.assert_array_almost_equal(
+        spl_jlkf.squeeze(),
+        [[2.21381744e+01, 2.60932044e+01, 2.88772888e+01, 2.84115953e+01, 2.82813953e+01, 2.55784584e+01,
+          7.29692738e+00, -5.38230264e+01],
+         [1.78459253e+01, 2.16729775e+01, 2.40686797e+01, 2.29153778e+01, 2.21888885e+01, 1.89631674e+01,
+          -3.52596583e-02, -6.18125090e+01]])
+    npt.assert_array_almost_equal(nmap['Total sound pressure level'].interp(x=[x[0], x[-1]], y=1000),
+                                  total_spl_jlk, 2)
+    npt.assert_array_almost_equal(nmap['Sound pressure level'].interp(x=[x[0], x[-1]], y=1000),
+                                  spl_jlkf, 1)
+
+    if 0:
+        import matplotlib.pyplot as plt
+        ax1, ax2, ax3 = plt.subplots(1, 3, figsize=(16, 4))[1]
+        sim_res.flow_map().plot_wake_map(ax=ax1)
+
+        ax1.plot([x[0]], [1000], '.', label='Receiver 1')
+        ax1.plot([x[-1]], [1000], '.', label='Receiver 2')
+        ax1.legend()
+
+        ax2.plot(nm.freqs, spl_jlkf[0, 0, 0], label='Receiver 1')
+        ax2.plot(nm.freqs, spl_jlkf[1, 0, 0], label='Receiver 2')
+        setup_plot(xlabel='Frequency [Hz]', ylabel='Sound pressure level [dB]', ax=ax2)
+        plt.tight_layout()
+
+        nmap['Total sound pressure level'].squeeze().plot(ax=ax3)
+        wt.plot(x, y, ax=ax3)
+
+        plt.show()

py_wake/wind_farm_models/wind_farm_model.py

@@ -12,6 +12,7 @@ import multiprocessing
 from py_wake.utils.parallelization import get_pool_map, get_pool_starmap
 from py_wake.utils.functions import arg2ilk, coords2ILK
 from py_wake.utils.gradients import autograd
+from py_wake.noise_models.iso import ISONoiseModel
 
 
 class WindFarmModel(ABC):
@@ -589,6 +590,25 @@ class SimulationResult(xr.Dataset):
 
         return ds
 
+    def noise_model(self, noiseModel=ISONoiseModel):
+        WS_eff_ilk = self.WS_eff_ilk
+        freqs, sound_power_level = self.windFarmModel.windTurbines.sound_power_level(WS_eff_ilk, **self.wt_inputs)
+        return noiseModel(src_x=self.x.values, src_y=self.y.values, src_h=self.h.values,
+                          freqs=freqs, sound_power_level=sound_power_level,
+                          elevation_function=self.windFarmModel.site.elevation)
+
+    def noise_map(self, noiseModel=ISONoiseModel, grid=None, ground_type=0, temperature=20, relative_humidity=80):
+        if grid is None:
+            grid = HorizontalGrid(h=2)
+        nm = self.noise_model(noiseModel)
+        X, Y, x_j, y_j, h_j, plane = self._get_grid(grid)
+        spl_jlk, spl_jlkf = nm(x_j, y_j, h_j, temperature, relative_humidity, ground_type=ground_type)
+        return xr.Dataset({'Total sound pressure level': (('y', 'x', 'wd', 'ws'),
+                                                          spl_jlk.reshape(X.shape + (spl_jlk.shape[1:]))),
+                           'Sound pressure level': (('y', 'x', 'wd', 'ws', 'freq'),
+                                                    spl_jlkf.reshape(X.shape + (spl_jlkf.shape[1:])))},
+                          coords={'x': X[0], 'y': Y[:, 0], 'wd': self.wd, 'ws': self.ws, 'freq': nm.freqs})
+
     def flow_box(self, x, y, h, wd=None, ws=None):
         X, Y, H = np.meshgrid(x, y, h)
         x_j, y_j, h_j = X.flatten(), Y.flatten(), H.flatten()

py_wake/wind_turbines/power_ct_functions.py

@@ -381,10 +381,11 @@ class PowerCtXr(PowerCtNDTabular):
             list of additional models.
         """
         assert method == 'linear'
-        assert 'power' in ds
-        assert 'ct' in ds
+        keys = {k.lower(): k for k in ds}
+        assert 'power' in keys
+        assert 'ct' in keys
         assert 'ws' in ds.dims
-        ds = ds[['power', 'ct']]
+        ds = ds[[keys['power'], keys['ct']]]
         power_arr, ct_arr = ds.to_array()
 
         if list(power_arr.dims).index('ws') > 0: