From 2b34b2fc7645bf48a6c0b0f5bb674eb4dccf752c Mon Sep 17 00:00:00 2001
From: Julian Quick <juqu@dtu.dk>
Date: Wed, 1 Jul 2026 22:43:24 -0700
Subject: [PATCH] fix(pywake): honor windIO flow-field bounds; test statistical
 flow field
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Flow-field output for statistical (wind rose) input works since the
run_pywake decomposition; add a regression test that runs a Weibull
rose and checks FarmFlow.nc contents.

The bounds override read z_planes keys xlb/xub/ylb/yub, which windIO
does not define — every example uses x_bounds/y_bounds, so requested
bounds were silently ignored and the flow box always fell back to the
site bounding box. Read the windIO keys instead.

Closes #7
---
 tests/test_pywake.py | 64 ++++++++++++++++++++++++++++++++++++++++++++
 wifa/pywake_api.py   |  9 +++----
 2 files changed, 68 insertions(+), 5 deletions(-)

diff --git a/tests/test_pywake.py b/tests/test_pywake.py
index 2605cf5..7e1dde5 100644
--- a/tests/test_pywake.py
+++ b/tests/test_pywake.py
@@ -466,6 +466,70 @@ def test_pywake_dict_timeseries_per_turbine_with_density(tmp_path):
     assert aep_with != aep_without
 
 
+def test_pywake_statistical_flow_field(tmp_path):
+    """Flow field output for statistical (wind rose) input (GH issue #7)."""
+    from conftest import _ANALYSIS, _TURBINE
+
+    system = {
+        "name": "statistical flow field test",
+        "site": {
+            "name": "Test site",
+            "boundaries": {
+                "polygons": [{"x": [-100, 800, 800, -100], "y": [100, 100, -100, -100]}]
+            },
+            "energy_resource": {
+                "name": "Test resource",
+                "wind_resource": {
+                    "wind_direction": [270.0],
+                    "wind_speed": list(range(4, 26)),
+                    "weibull_a": {"data": [9.0], "dims": ["wind_direction"]},
+                    "weibull_k": {"data": [2.2], "dims": ["wind_direction"]},
+                    "sector_probability": {"data": [1.0], "dims": ["wind_direction"]},
+                    "turbulence_intensity": {
+                        "data": [0.07],
+                        "dims": ["wind_direction"],
+                    },
+                },
+            },
+        },
+        "wind_farm": {
+            "name": "Test farm",
+            "layouts": [{"coordinates": {"x": [0.0, 700.0], "y": [0.0, 0.0]}}],
+            "turbines": _TURBINE,
+        },
+        "attributes": {
+            "flow_model": {"name": "pywake"},
+            "analysis": _ANALYSIS,
+            "model_outputs_specification": {
+                "flow_field": {
+                    "report": True,
+                    "output_variables": ["velocity_u", "turbulence_intensity"],
+                    "z_planes": {
+                        "z_sampling": "hub_heights",
+                        "xy_sampling": "grid",
+                        "x_bounds": [-500.0, 1500.0],
+                        "y_bounds": [-400.0, 400.0],
+                        "dx": 100.0,
+                        "dy": 100.0,
+                    },
+                },
+            },
+        },
+    }
+
+    run_pywake(system, output_dir=str(tmp_path))
+
+    flow_file = tmp_path / "FarmFlow.nc"
+    assert flow_file.exists()
+    with xr.open_dataset(flow_file) as ds:
+        assert {"wind_speed", "turbulence_intensity"} <= set(ds.data_vars)
+        # requested x/y bounds are honored (windIO x_bounds/y_bounds keys)
+        assert float(ds.x.min()) == -500.0 and float(ds.x.max()) >= 1400.0
+        assert float(ds.y.min()) == -400.0 and float(ds.y.max()) >= 300.0
+        # flow field sampled at hub height
+        assert list(ds.z.values) == [_TURBINE["hub_height"]]
+
+
 # if __name__ == "__main__":
 #    test_heterogeneous_wind_rose()
 #     simple_yaml_to_pywake('../examples/cases/windio_4turbines_multipleTurbines/plant_energy_turbine/IEA_10MW_turbine.yaml')
diff --git a/wifa/pywake_api.py b/wifa/pywake_api.py
index 91cf99a..ffa05bb 100644
--- a/wifa/pywake_api.py
+++ b/wifa/pywake_api.py
@@ -240,7 +240,7 @@ def get_flow_field_param(system_dat, param_name, default=None):
 
     Args:
         system_dat: System data dictionary
-        param_name: Name of parameter to extract (e.g., 'xlb', 'dx')
+        param_name: Name of parameter to extract (e.g., 'x_bounds', 'dx')
         default: Default value if parameter not found
 
     Returns:
@@ -279,10 +279,9 @@ def construct_site(system_dat, resource_dat, hub_heights, x_positions):
     WFYUB = np.max(boundaries["y"])
 
     # Override with explicit flow field bounds if specified
-    WFXLB = get_flow_field_param(system_dat, "xlb", WFXLB)
-    WFXUB = get_flow_field_param(system_dat, "xub", WFXUB)
-    WFYLB = get_flow_field_param(system_dat, "ylb", WFYLB)
-    WFYUB = get_flow_field_param(system_dat, "yub", WFYUB)
+    # (windIO z_planes uses x_bounds/y_bounds pairs)
+    WFXLB, WFXUB = get_flow_field_param(system_dat, "x_bounds", (WFXLB, WFXUB))
+    WFYLB, WFYUB = get_flow_field_param(system_dat, "y_bounds", (WFYLB, WFYUB))
     WFDX = get_flow_field_param(system_dat, "dx", (WFXUB - WFXLB) / 100)
     WFDY = get_flow_field_param(system_dat, "dy", (WFYUB - WFYLB) / 100)