diff --git a/examples/examples_floating/002_floating_vs_fixedbottom_farm.py b/examples/examples_floating/002_floating_vs_fixedbottom_farm.py
index 0400ac7f1d..bfcd683b59 100644
--- a/examples/examples_floating/002_floating_vs_fixedbottom_farm.py
+++ b/examples/examples_floating/002_floating_vs_fixedbottom_farm.py
@@ -24,8 +24,6 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
-import pandas as pd
-from scipy.interpolate import NearestNDInterpolator
 
 import floris.flow_visualization as flowviz
 from floris import FlorisModel, WindRose
diff --git a/floris/core/core.py b/floris/core/core.py
index 1b78423f2a..d64cc91c45 100644
--- a/floris/core/core.py
+++ b/floris/core/core.py
@@ -283,6 +283,11 @@ def solve_for_velocity_deficit_profiles(
         for more details.
         """
 
+        self.logger.warning(
+            "Velocity deficit profiles will move to a Numpy data structure in the next release. "
+            "See https://github.com/NatLabRockies/floris/pull/1194."
+        )
+
         # Create a grid that contains coordinates for all the sample points in all profiles.
         # Effectively, this is a grid of parallel lines.
         n_lines = len(downstream_dists)
diff --git a/floris/core/turbine/turbine.py b/floris/core/turbine/turbine.py
index d9900037ea..900feca913 100644
--- a/floris/core/turbine/turbine.py
+++ b/floris/core/turbine/turbine.py
@@ -6,7 +6,6 @@
 
 import attrs
 import numpy as np
-import pandas as pd
 from attrs import define, field
 from scipy.interpolate import interp1d
 
@@ -656,34 +655,40 @@ def _initialize_multidim_power_thrust_table(self):
         self.power_thrust_data_file = self.turbine_library_path / self.power_thrust_data_file
 
         # Read in the multi-dimensional data supplied by the user.
-        df = pd.read_csv(self.power_thrust_data_file)
+        data = np.genfromtxt(self.power_thrust_data_file, delimiter=',', names=True)
 
-        # Down-select the DataFrame to have just the ws, Cp, and Ct values
-        index_col = df.columns.values[:-3]
-        self.condition_keys = index_col.tolist()
-        df2 = df.set_index(index_col.tolist())
+        # The CSV columns are: [condition_keys..., ws, power, thrust_coefficient]
+        # Condition keys are the leading columns (e.g. Tp, Hs)
+        # ws/power/thrust_coefficient are always the last 3 and are excluded.
+        self.condition_keys = list(data.dtype.names[:-3])
+
+        # Find unique combinations of condition key values.
+        # np.column_stack promotes 1D arrays to (N, 1), so this works for any number of keys.
+        cond_data = np.column_stack([data[c] for c in self.condition_keys])
+        unique_keys = [tuple(row) for row in np.unique(cond_data, axis=0)]
 
         # Loop over the multi-dimensional keys to get the correct ws/Cp/Ct data to make
         # the thrust_coefficient and power interpolants.
         power_thrust_table_ = {} # Reset
-        for key in df2.index.unique():
-            # Select the correct ws/Cp/Ct data
-            data = df2.loc[key]
-            if type(key) is not tuple:
-                key = (key,)
+        for key in unique_keys:
+            # Build a boolean mask selecting rows that match this condition combination
+            mask = np.ones(len(data), dtype=bool)
+            for col, val in zip(self.condition_keys, key):
+                mask &= data[col] == val
+
+            rows = data[mask]
 
             # Build the interpolants
             power_thrust_table_.update(
                 {
                     key: {
-                        "wind_speed": data['ws'].values,
-                        "power": data['power'].values,
-                        "thrust_coefficient": data['thrust_coefficient'].values,
+                        "wind_speed": rows['ws'],
+                        "power": rows['power'],
+                        "thrust_coefficient": rows['thrust_coefficient'],
                         **power_thrust_table_ref
                     },
                 }
             )
-            # Add reference information at the lower level
 
         # Save names of dimensions and set on-object version
         power_thrust_table_.update({"condition_keys": self.condition_keys})
diff --git a/pyproject.toml b/pyproject.toml
index 30bc376d25..d6860da05a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -36,7 +36,7 @@ dependencies = [
     "numpy~=2.0",
     "scipy~=1.1",
     "matplotlib~=3.0",
-    "pandas~=2.0",
+    "pandas>=2.0,<4",
     "shapely~=2.0",
     "coloredlogs~=15.0",
     "pathos~=0.3",
diff --git a/tests/geometric_yaw_unit_test.py b/tests/geometric_yaw_unit_test.py
index 61edafc45d..e32d2be1b1 100644
--- a/tests/geometric_yaw_unit_test.py
+++ b/tests/geometric_yaw_unit_test.py
@@ -1,6 +1,5 @@
 
 import numpy as np
-import pandas as pd
 
 from floris import FlorisModel
 from floris.optimization.yaw_optimization.yaw_optimizer_geometric import (
diff --git a/tests/reg_tests/random_search_layout_opt_regression_test.py b/tests/reg_tests/random_search_layout_opt_regression_test.py
index d5f3313de2..3725b99273 100644
--- a/tests/reg_tests/random_search_layout_opt_regression_test.py
+++ b/tests/reg_tests/random_search_layout_opt_regression_test.py
@@ -1,6 +1,5 @@
 
 import numpy as np
-import pandas as pd
 
 from floris import FlorisModel, WindRose
 from floris.optimization.layout_optimization.layout_optimization_random_search import (
diff --git a/tests/reg_tests/scipy_layout_opt_regression.py b/tests/reg_tests/scipy_layout_opt_regression.py
index 1029dfd764..2d10a64ac6 100644
--- a/tests/reg_tests/scipy_layout_opt_regression.py
+++ b/tests/reg_tests/scipy_layout_opt_regression.py
@@ -1,6 +1,5 @@
 
 import numpy as np
-import pandas as pd
 
 from floris import FlorisModel, WindRose
 from floris.optimization.layout_optimization.layout_optimization_scipy import (
diff --git a/tests/serial_refine_unit_test.py b/tests/serial_refine_unit_test.py
index cfda030a7c..98deca2c3f 100644
--- a/tests/serial_refine_unit_test.py
+++ b/tests/serial_refine_unit_test.py
@@ -1,6 +1,5 @@
 
 import numpy as np
-import pandas as pd
 
 from floris import FlorisModel
 from floris.optimization.yaw_optimization.yaw_optimizer_sr import YawOptimizationSR
diff --git a/tests/turbine_multi_dim_unit_test.py b/tests/turbine_multi_dim_unit_test.py
index 0a0f2882a2..11e4ea10b4 100644
--- a/tests/turbine_multi_dim_unit_test.py
+++ b/tests/turbine_multi_dim_unit_test.py
@@ -2,7 +2,6 @@
 from pathlib import Path
 
 import numpy as np
-import pandas as pd
 import pytest
 
 from floris.core import (
@@ -23,29 +22,6 @@
 
 INDEX_FILTER = [0, 2]
 
-# NOTE: MultiDimensionalPowerThrustTable not used anywhere, so I'm commenting
-# this out.
-
-# def test_multi_dimensional_power_thrust_table():
-#     turbine_data = SampleInputs().turbine_multi_dim
-#     turbine_data["power_thrust_data_file"] = CSV_INPUT
-#     df_data = pd.read_csv(turbine_data["power_thrust_data_file"])
-#     flattened_dict = MultiDimensionalPowerThrustTable.from_dataframe(df_data)
-#     flattened_dict_base = {
-#         ('Tp', '2', 'Hs', '1'): [],
-#         ('Tp', '2', 'Hs', '5'): [],
-#         ('Tp', '4', 'Hs', '1'): [],
-#         ('Tp', '4', 'Hs', '5'): [],
-#     }
-#     assert flattened_dict == flattened_dict_base
-
-#     # Test for initialization errors
-#     for el in ("ws", "Cp", "Ct"):
-#         df_data = pd.read_csv(turbine_data["power_thrust_data_file"])
-#         df = df_data.drop(el, axis=1)
-#         with pytest.raises(ValueError):
-#             MultiDimensionalPowerThrustTable.from_dataframe(df)
-
 
 def test_turbine_init():
     turbine_data = SampleInputs().turbine_multi_dim
diff --git a/tests/yaw_optimization_integration_test.py b/tests/yaw_optimization_integration_test.py
index a0c3011fc8..7bbcc23dda 100644
--- a/tests/yaw_optimization_integration_test.py
+++ b/tests/yaw_optimization_integration_test.py
@@ -1,5 +1,4 @@
 import numpy as np
-import pandas as pd
 import pytest
 
 from floris import FlorisModel