fi.calculate_horizontal_plane wind_direction angle problem

[acse-xm421]

Here is part of my code:

  # Create the plots of un-optimized layout
  show_wd = 30
  horizontal_plane = fi.calculate_horizontal_plane(
      x_resolution=200,
      y_resolution=100,
      height=90.0,
      yaw_angles=yaw_angles,
      wd=[show_wd],
      ws=[8.0],
  )
  fig, ax_list = plt.subplots(1, 1, figsize=(10, 8))
  wakeviz.visualize_cut_plane(horizontal_plane, ax=ax_list, title="Horizontal")
  name = "examples/test_pic/" + save_name + "_" + str(show_wd) + ".png"
  plt.savefig(name)

But the plot looks like this.

Only when show_wd = 0 or 90 or something else, it works well.

[rafmudaf]

Hi @acse-xm421 can you describe what you expected to see in the plot with this wind direction?

[RHammond2]

I think I might have a simple fix to this that could be implemented in the plotting routine by extending the plane to the corners, and using the u value of the nearest point and plotting the countour from there. It might get weird for heterogeneous, but if we really wanted to, I think we could extend the corners to a series of points on the frame's edge.

[RHammond2]

@acse-xm421 the below should be able to be incorporated into your code as needed without any additional package installations.

@rafmudaf / @paulf81 let me know if the below is an additional we'd like to see in the code with calculate_fill and use_ambient kwargs and a check for if the corners already exist in the plane (proxy for rotated plane). I'm sure there's also a better way to grab the ambient wind speed, but I'm not totally sure what that is without really digging deep.

# Modififed example 2

import matplotlib.pyplot as plt
import numpy as np

import floris.tools.visualization as wakeviz
from floris.tools import FlorisInterface
from floris.tools.visualization import (
    calculate_horizontal_plane_with_turbines,
    visualize_cut_plane,
)

# NOTE: New imports below
from copy import deepcopy

import pandas as pd
from scipy import spatial, stats

# Get the FI object and calculate the horizontal plane, rotated at 30 degrees (should work for any points)
fi = FlorisInterface("inputs/gch.yaml")
horizontal_plane = fi.calculate_horizontal_plane(
    x_resolution=200,
    y_resolution=100,
    height=90.0,
    yaw_angles=np.array([[[25.,0.,0.]]]),
    wd=[30],
)

# Calculate the new plane
use_ambient = False
plane_w_frame = horizontal_plane.df.copy()

# 1. Retrieve the coordinate pairs and find the corners of the grid
coords = np.array(list(zip(horizontal_plane.df.x1, horizontal_plane.df.x2)))
lower_left = coords.min(axis=0)
upper_right = coords.max(axis=0)
lower_right = np.array([upper_right[0], lower_left[1]])
upper_left = np.array([lower_left[0], upper_right[1]])


# 2. If one corner is included, they're all included due to the nature of rotational grids
# NOTE: a check for the inclusion of a corner goes here, but isn't needed for this example
tree = spatial.KDTree(plane_w_frame[["x1", "x2"]])

# 3. Create points along the frame edge (using 101 points to include endpoints and have sufficient coverage)
right_frame_coords = np.hstack((np.full(101, upper_right[0]).reshape(-1, 1), np.linspace(upper_right[1], lower_right[1], 101).reshape(-1, 1)))
left_frame_coords = np.hstack((np.full(101, upper_left[0]).reshape(-1, 1), np.linspace(upper_left[1], lower_left[1], 101).reshape(-1, 1)))
top_frame_coords = np.hstack((np.linspace(upper_left[0], upper_right[0], 101).reshape(-1, 1), np.full(101, upper_right[1]).reshape(-1, 1)))
lower_frame_coords = np.hstack((np.linspace(lower_left[0], lower_right[0], 101).reshape(-1, 1), np.full(101, lower_right[1]).reshape(-1, 1)))
frame_coords = np.vstack((top_frame_coords, right_frame_coords, lower_frame_coords, left_frame_coords))

# 4. Get the closest value to those points and repopulate the cut plane
frame_nn_ix = tree.query(frame_coords)[1]
u = plane_w_frame.iloc[frame_nn_ix].u.values
u = u if not use_ambient else np.full(u.shape, stats.mode(u, keepdims=False).mode)
frame_df = pd.DataFrame(frame_coords, columns=["x1", "x2"])
frame_df = frame_df.assign(
    x3=plane_w_frame.iloc[frame_nn_ix].x3.values,
    u=u,
    v=plane_w_frame.iloc[frame_nn_ix].v.values,
    w=plane_w_frame.iloc[frame_nn_ix].w.values,
)
plane_w_frame = pd.concat([plane_w_frame, frame_df], axis=0).reset_index(drop=True)

# 5. Plot the new results with our updated horizontal plane
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111)

modified_horizontal_plane = deepcopy(horizontal_plane)
modified_horizontal_planedf = plane_w_frame

wakeviz.visualize_cut_plane(modified_horizontal_plane, ax=ax, title="Horizontal")
fig.show()

Results for use_ambient=True

Results for use_ambient=False. Notice here how the bounds in the original below cut off the end of the wake field, so there is some weirdness where those are the nearest points to those on the frame.

Original

[rafmudaf]

For the second case could you try extending the bounds downstream? The default is 10D so if you increase that to 15D or so it might resolve those numerical problems. I don't remember when the rotation happens, so I would try first to extend y_bounds and then if that doesn't work you may have to extend both y_bounds and x_bounds.

[rafmudaf]

Otherwise, I like the result here, so good work @RHammond2. We should think through how this fits in to the rest of the visualization tools, and also how it works with heterogeneous wind conditions. As I understand it, this is a post processing step for visualization only, but I wonder if we shouldn't extend the flow field that is calculated in floris and then clip it to fit into a matplotlib window. It would be a similar result but from a different method. The one above seems simpler and more straightforward, but the approach I mention here may get around any problems with the heterogeneous inputs.

This general topic was discussed during #578. If I remember right, the decision there was to set the areas outside the bounds of the heterogeneous wind condition definition to some very low value. Given that, maybe the decision here should be that this is only valid for a constant wind condition.

[RHammond2]

Thanks for that @rafmudaf! Seems like I need to get to know the ins and outs of the tool library a bit better. So this technically works better for the results, but is harder to manage from a plotting perspective, and would require some consideration on the users' part.

The below, without any of the fancy stuff from the previous example will produce nearly identical results, but notice that there is now guesswork in producing a plot of the area neatly, and that if we don't set the axis limits, we end up with the original result, but a bigger bound. Ultimately, I think it's the better solution, but in a future tools/visualization refactor, we could consider allowing user provided axis limits. Ideally I think we'd have a way to plot this more consistently, but that will be quite tricky, and it might be better to show users how to get the results they want in a short plotting tutorial, rather than coming up with a complex solution for finding the idealized plotting window.

fi = FlorisInterface("inputs/gch.yaml")
buffer = 15 * fi.floris.farm.rotor_diameters.mean()
xmin = fi.floris.farm.layout_x.min() - 2 * buffer
xmax = fi.floris.farm.layout_x.max() + 2 * buffer
ymin = fi.floris.farm.layout_y.min() - 2 * buffer
ymax = fi.floris.farm.layout_y.max() + 2 * buffer
horizontal_plane = fi.calculate_horizontal_plane(
    x_resolution=200,
    y_resolution=100,
    x_bounds=(xmin, xmax),
    y_bounds=(ymin, ymax),
    height=90.0,
    yaw_angles=np.array([[[25.,0.,0.]]]),
    wd=[30],
)

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111)

wakeviz.visualize_cut_plane(horizontal_plane, ax=ax, title="Horizontal")
ax.set_xlim(xmin + buffer, xmax - 1.5 * buffer)
ax.set_ylim(ymin + 0.5 * buffer, ymax - 1.5 * buffer)
fig.show()

With custom axis limits

Without axis limits