Coverage for wifa_uq / model_error_database / utils.py: 56%
112 statements
« prev ^ index » next coverage.py v7.13.0, created at 2025-12-19 02:10 +0000
1import numpy as np
2import matplotlib.pyplot as plt
3from scipy.spatial import ConvexHull
6def calc_boundary_area(x, y, show=True):
7 """
8 Use convex hull to calculate area of wind farm boundary and optionally plot
9 """
10 points = np.column_stack((x, y))
11 hull = ConvexHull(points)
13 if show:
14 plt.scatter(x, y, color="blue", label="Turbines")
15 # Plot boundary
16 for simplex in hull.simplices:
17 plt.plot(points[simplex, 0], points[simplex, 1], "r-")
19 # Close the hull
20 plt.plot(
21 [points[hull.vertices[0], 0], points[hull.vertices[-1], 0]],
22 [points[hull.vertices[0], 1], points[hull.vertices[-1], 1]],
23 "r-",
24 )
25 plt.xlabel("x-coordinate")
26 plt.ylabel("y-coordinate")
27 plt.title(f"Wind Farm Boundary (Area = {hull.volume:.2f})")
28 plt.legend()
29 plt.axis("equal")
30 plt.show()
32 return hull.volume
35def blockage_metrics(
36 xy, wind_dir_deg_from_north, D, grid_res=151, L_inf_factor=20.0, plot=False
37):
38 """
39 Compute per-turbine Blockage Ratio (BR) and Blockage Distance (BD)
40 for a given 2D wind farm layout and wind direction.
42 Parameters
43 ----------
44 xy : (N, 2) array
45 Turbine coordinates in meters (X east, Y north).
46 wind_dir_deg_from_north : float
47 Meteorological wind direction in degrees (0 = from north, 90 = from east).
48 D : float
49 Rotor diameter (meters). All turbines assumed identical (as in the paper).
50 grid_res : int
51 Number of samples per axis for disk discretization (odd number recommended).
52 Effective number of points on the rotor ≈ π*(grid_res/2)^2 / 4.
53 L_inf_factor : float
54 L∞ as a multiple of D (paper uses 20D).
56 Returns
57 -------
58 BR : (N,) array
59 Blockage ratio for each turbine (0..1).
60 BD : (N,) array
61 Blockage distance for each turbine (meters).
62 BR_farm : float
63 Mean BR across turbines.
64 BD_farm : float
65 Mean BD across turbines.
66 """
68 xy = np.asarray(xy, dtype=float)
69 N = xy.shape[0]
70 R = 0.5 * D
72 # Default blockage distance for turbines which aren't blocked
73 L_inf = L_inf_factor * D
75 # Change direction into radians and define downwind and crosswind directions
76 theta = np.deg2rad(wind_dir_deg_from_north)
77 a_hat = np.array([-np.sin(theta), -np.cos(theta)])
78 n_hat = np.array([np.cos(theta), -np.sin(theta)])
80 # map turbine coordinates to downwind and crosswind vectors
81 x_down = xy @ a_hat
82 x_cross = xy @ n_hat
84 # Creating a mesh grid for a rotor disk in the (u,z) plane:
85 # u=horizontal crosswind offset, z=vertical offset (offsets from the centre of disk)
86 # Keeping only points inside the rotor disk
87 lin = np.linspace(-R, R, grid_res)
88 uu, zz = np.meshgrid(lin, lin, indexing="xy")
89 mask = (uu * uu + zz * zz) <= (R * R)
90 u_pts = uu[mask]
91 z_pts = zz[mask]
92 M = u_pts.size
94 # Initializing arrays to store blockage ratio and blockage distance for each turbine
95 BR = np.zeros(N)
96 BD = np.zeros(N)
98 # For all turbines
99 for i in range(N):
100 # Initialize the grid for our turbine disk
101 u_i = u_pts
102 # z_i = z_pts
104 # calculate the downstream distance between this turbine and all others
105 delta_down = x_down[i] - x_down # (N,)
106 upstream_mask = delta_down > 0.0
107 # if a turbine is not blocked
108 if not np.any(upstream_mask):
109 BR[i] = 0.0
110 BD[i] = 1.0
111 continue
113 # For upstream turbines, calculate the crosswind and downwind distance between rotor centres
114 cross_dist = x_cross[i] - x_cross[upstream_mask]
115 down_dist = delta_down[upstream_mask]
117 ### Initialize some variables:
118 # boolean array representing which grid cells of the rotor grid are blocked by any upstream turbine
119 blocked = np.zeros(M, dtype=bool)
121 # array with default value of L_inf containing the distance to the nearest upstream turbine
122 L_point = np.full(M, L_inf, dtype=float)
124 ### Looping over each upstream turbine and calculating blockage on all rotor grid cells
125 for j in range(cross_dist.size):
126 dx = u_i - cross_dist[j] # horizontal distance to upstream turbine centre
127 dz = z_pts # vertical offset
128 d = np.sqrt(dx**2 + dz**2) # actual straight line /euclidean distance
129 hits = (
130 d <= R
131 ) # boolean array of rotor grid cells that are blocked by this upstream turbine
133 if np.any(hits):
134 blocked[hits] = True
135 # In the scenario that a grid point is blocked by multiple turbines, we select the nearest one
136 # This is done by selecting the minimum distance from either the turbine in the previous iteration, or the
137 # Distance to the turbine in the current iteration
138 L_point[hits] = np.minimum(L_point[hits], down_dist[j])
140 # Calculating the ratios over the disk
141 BR[i] = blocked.mean()
143 # "blocked" is a boolean representing which grid cells are blocked and L_point represents the distance from each blocked point to the nearest blocking turbine
144 # If there is no grid cells blocked, we default to L_inf and the total metric is a weighted sum dependent on how many grid cells are blocked
145 BD[i] = (blocked * L_point + (~blocked) * L_inf).mean() / L_inf
147 # Quantities on a farm level
148 BR_farm = BR.mean()
149 BD_farm = BD.mean()
151 if plot:
153 def plot_metric(metric, title, cmap="viridis"):
154 fig, ax = plt.subplots(figsize=(8, 8))
155 sc = ax.scatter(
156 xy[:, 0], xy[:, 1], c=metric, s=100, cmap=cmap, edgecolor="k"
157 )
158 cbar = plt.colorbar(sc, ax=ax)
159 cbar.set_label("Metric value", fontsize=12)
161 # Wind arrow: points FROM where wind is coming
162 arrow_length = max(np.ptp(xy[:, 0]), np.ptp(xy[:, 1])) * 0.1
163 ax.arrow(
164 np.mean(xy[:, 0]),
165 np.mean(xy[:, 1]),
166 -arrow_length * np.sin(theta),
167 -arrow_length * np.cos(theta),
168 head_width=arrow_length * 0.2,
169 head_length=arrow_length * 0.2,
170 fc="r",
171 ec="r",
172 linewidth=2,
173 )
174 ax.text(
175 np.mean(xy[:, 0]),
176 np.mean(xy[:, 1]) + arrow_length * 0.15,
177 "Wind",
178 color="r",
179 ha="center",
180 )
182 pad = max(250, 0.1 * max(np.ptp(xy[:, 0]), np.ptp(xy[:, 1])))
183 x_min, x_max = xy[:, 0].min() - pad, xy[:, 0].max() + pad
184 y_min, y_max = xy[:, 1].min() - pad, xy[:, 1].max() + pad
185 ax.set_xlim(x_min, x_max)
186 ax.set_ylim(y_min, y_max)
188 ax.set_xlabel("X [m]")
189 ax.set_ylabel("Y [m]")
190 ax.set_aspect("equal")
191 ax.set_title(title)
192 plt.show()
194 plot_metric(
195 BR,
196 f"Blockage Ratio per Turbine: Farm Average = {BR_farm:.2f}",
197 cmap="viridis",
198 )
199 plot_metric(
200 BD,
201 f"Blockage Distance per Turbine [m]: Farm Average = {BD_farm:.2f}",
202 cmap="plasma",
203 )
205 return BR, BD, BR_farm, BD_farm
208def farm_length_width(x, y, wind_dir_deg_from_north, D, plot=False):
209 """
210 Compute farm length and width in wind and crosswind directions, normalized by diameter.
212 """
213 xy = np.column_stack((x, y))
214 # N = len(xy)
215 theta = np.deg2rad(wind_dir_deg_from_north)
216 # Wind direction unit vector
217 wind_vec = np.array([-np.sin(theta), -np.cos(theta)])
218 cross_vec = np.array([np.cos(theta), -np.sin(theta)])
220 # Project all points onto wind and crosswind axes
221 proj_wind = xy @ wind_vec
222 proj_cross = xy @ cross_vec
224 length = np.round((proj_wind.max() - proj_wind.min()) / D)
225 width = np.round((proj_cross.max() - proj_cross.min()) / D)
227 # Optional plotting
228 if plot:
230 def plot_turbine_layout(x, y, wind_dir_deg_from_north, title="Turbine Layout"):
231 plt.figure(figsize=(6, 6))
232 plt.scatter(x, y, s=100, c="b", label="Turbines")
233 plt.xlabel("x [m]")
234 plt.ylabel("y [m]")
235 plt.axis("equal")
236 plt.grid(True)
237 plt.title(title)
239 # Wind direction arrow
240 center_x = np.mean(x)
241 center_y = np.mean(y)
242 arrow_length = max(np.ptp(x), np.ptp(y)) * 0.2
243 theta = np.deg2rad(wind_dir_deg_from_north)
244 dx = -arrow_length * np.sin(theta)
245 dy = -arrow_length * np.cos(theta)
246 plt.arrow(
247 center_x,
248 center_y,
249 dx,
250 dy,
251 head_width=arrow_length * 0.15,
252 head_length=arrow_length * 0.15,
253 fc="r",
254 ec="r",
255 linewidth=2,
256 )
257 plt.text(
258 center_x + dx * 1.1,
259 center_y + dy * 1.1,
260 "Wind",
261 color="r",
262 ha="center",
263 va="center",
264 )
265 plt.legend()
266 plt.show()
268 plot_turbine_layout(
269 x,
270 y,
271 wind_dir_deg_from_north=wind_dir_deg_from_north,
272 title=f"D= {D}\n Length={length}D\nWidth={width}D",
273 )
275 return length, width