Add WRG file support to FLORIS (#919)

* Add wrg_reader file * Blocking in notebook to generate example wrg file * small additions * Update notebook and example wrg * Remove wrg reader file * Notebook to script * add WindResourceGrid to import * Initial implementation of WindResourceGrid * First tests of wind_resource_grid * Change to 2x3 * Update wind_resource_grid and tests * Add print out * Add print out * Add WindRoseByTurbineObject * Method to generate WindRoseByTurbine * Add tests of WindRoseByTurbine * test get_wind_rose_by_turbine * Update examples * Add a test that expected farm power hasn't changed * Start implementing expected turbine power * Update example * Update unpack * bugfix * Update expected functions * Finish example 2 * Add test for expected functions * Renumber examples * Add wrg option to layout opt * Add layout opt example * Update 001 to make clearer accelerations * Add print out feature * Add get_hetergeneous_map func * Add examples 4 and 5 * Update test * Remove echo from example * Renumber examples * Update init to new structure * Update example 001 * Delete wind_resource_grid file * Update to WindRoseWRG * Update imports and references to old class name * Update example 1 * Update example 01 * Update example 001 * Update example 1 * Floris updates layout * Fix update layout * Allow ti_table * fix call to set_layout * Update example 2 * Fix the ordering of applying layout to wind data * Roll back changes to layout opt * Update example 3 * clean up plot * back to 60 * Fix issue in upsampling limits * Update example 1 * Rename example 2 * Update example 3 * Bugfix in random search * Update layout examples * Update example wrg * Update het function * Update tests * Update example 3 * A couple of minor formatting changes. * Improve upsampling to retain coverage of original areas * Update tests * Initial upsample fixes * Upsample testing changes * clean up code * clean up tests * matching changes to wind_ti_rose * update wind_ti_rose tests * Split off get_weighted_turbine_powers function * Update test names * Fix call to equal axis * Update wind_data doc * Small update * Formatting updates. * Minor cleanup. * Remove subset_wind_speeds from WindRose and perform actions in WindRoseWRG.get_heterogeneous_wind_rose instead. * bugfix * Remove block comment * Add and test check for even spacing of wind_speeds; add note about assigned freq. --------- Co-authored-by: misi9170 <michael.sinner@nrel.gov>
NREL · Jul 30, 2024 · bb9e1b5 · bb9e1b5
1 parent 296628a
commit bb9e1b5
Show file tree

Hide file tree

Showing 12 changed files with 2,055 additions and 151 deletions.
diff --git a/docs/wind_data_user.ipynb b/docs/wind_data_user.ipynb
diff --git a/examples/examples_wind_resource_grid/000_generate_example_wrg.py b/examples/examples_wind_resource_grid/000_generate_example_wrg.py
@@ -0,0 +1,227 @@
+"""Example: Generate Example WRG File
+
+This first example demonstrates the content and structure of a
+Wind Resource Grid (WRG) file.
+
+WRG files are Wind Resource Grid files, and their structure is
+defined here:
+https://backend.orbit.dtu.dk/ws/portalfiles/portal/116352660/WAsP_10_Help_Facility.pdf
+
+In the script, a synthetic WRG file is derived using the WindRose class.
+
+"""
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.optimize import curve_fit
+
+from floris import WindRose
+
+
+# Define a function given the distribution of wind speeds in one sector,
+#  compute the A and k parameters of the Weibull distribution
+def weibull_func(U, A, k):
+    return (k / A) * (U / A) ** (k - 1) * np.exp(-((U / A) ** k))
+
+
+def estimate_weibull(U, freq):
+    # Normalize the frequency
+    freq = freq / freq.sum()
+
+    # Fit the Weibull distribution
+    popt, _ = curve_fit(weibull_func, U, freq, p0=(6.0, 2.0))
+    A_fit, k_fit = popt
+
+    return A_fit, k_fit
+
+
+##################################################
+# Parameters
+##################################################
+# Top line parameters
+Nx = 2  # Number of grid points in x
+Ny = 3  # Number of grid points in y
+Xmin = 0.0  # Minimum value of x (m)
+Ymin = 0.0  # Minimum value of y (m)
+cell_size = 1000.0  # Grid spacing (m)
+
+# Other fixed parameters
+z_coord = 0.0  # z-coordinate of the grid
+height_above_ground_level = 90.0  # Height above ground level
+num_sectors = 12  # Number of direction sectors
+
+
+
+##################################################
+# Generating data
+##################################################
+# The above parameters define a 3x3 grid of points.  Let's start
+# by assuming the point at (0,0) has the wind rose as
+# defined in inputs/wind_rose.csv
+wind_rose_base = WindRose.read_csv_long(
+    "../inputs/wind_rose.csv", wd_col="wd", ws_col="ws", freq_col="freq_val", ti_col_or_value=0.06
+)
+
+# Resample to number of sectors
+wind_rose_base = wind_rose_base.aggregate(wd_step=360 / num_sectors)
+
+## Generate the other wind roses
+# Assume that the wind roses at other points are generated by increasing
+# the north winds with increasing y and east winds with increasing x
+x_list = []
+y_list = []
+wind_rose_list = []
+
+for xi in range(Nx):
+    for yi in range(Ny):
+        # Get the x and y locations for this point
+        x = Xmin + xi * cell_size
+        y = Ymin + yi * cell_size
+        x_list.append(x)
+        y_list.append(y)
+
+        # Instantiate the wind rose object
+        wind_rose = WindRose.read_csv_long(
+            "../inputs/wind_rose.csv",
+            wd_col="wd",
+            ws_col="ws",
+            freq_col="freq_val",
+            ti_col_or_value=0.06,
+        )
+
+        # Resample to number of sectors
+        wind_rose = wind_rose.aggregate(wd_step=360 / num_sectors)
+
+        # Copy the frequency table
+        freq_table = wind_rose.freq_table.copy()
+
+        # How much to shift the wind rose for this location
+        percent_x = xi / (Nx - 1)
+        percent_y = yi / (Ny - 1)
+
+        # East frequency scaling
+        east_row = freq_table[3, :]
+        shift_amount = percent_x * east_row[:5] * .9
+        east_row[:5] = east_row[:5] - shift_amount
+        east_row[5:10] = east_row[5:10] + shift_amount
+        freq_table[3, :] = east_row
+
+        # North frequency scaling
+        north_row = freq_table[0, :]
+        shift_amount = percent_y * north_row[:6] * .9
+        north_row[:6] = north_row[:6] - shift_amount
+        north_row[6:12] = north_row[6:12] + shift_amount
+        freq_table[0, :] = north_row
+
+        # Add to list
+        wind_rose_list.append(
+            WindRose(
+                wind_directions=wind_rose.wind_directions,
+                wind_speeds=wind_rose.wind_speeds,
+                ti_table=wind_rose.ti_table,
+                freq_table=freq_table,
+            )
+        )
+
+##################################################
+# Show the wind roses in a grid
+##################################################
+
+fig, axarr = plt.subplots(Nx, Ny, figsize=(12, 12), subplot_kw={"polar": True})
+axarr = axarr.flatten()
+
+for i, wind_rose in enumerate(wind_rose_list):
+    wind_rose.plot(ax=axarr[i], ws_step=5)
+    axarr[i].set_title(f"({x_list[i]}, {y_list[i]})")
+
+fig.suptitle("Wind Roses at Grid Points")
+
+
+##################################################
+# Demonstrate fitting the Weibull distribution
+##################################################
+
+freq_test = wind_rose_list[0].freq_table[0, :] / wind_rose_list[0].freq_table[0, :].sum()
+a_test, k_test = estimate_weibull(wind_rose_list[0].wind_speeds, freq_test)
+print(f"A: {a_test}, k: {k_test}")
+
+fig, ax = plt.subplots(1, 1, figsize=(6, 3))
+ax.plot(wind_rose_list[0].wind_speeds, freq_test, label="Original")
+ax.plot(
+    wind_rose_list[0].wind_speeds,
+    weibull_func(wind_rose_list[0].wind_speeds, a_test, k_test),
+    label="Fitted",
+)
+ax.legend()
+ax.set_xlabel("Wind speed (m/s)")
+ax.set_ylabel("Frequency")
+ax.grid(True)
+
+
+##################################################
+# Write out the WRG file
+##################################################
+
+# Open the file
+with open("wrg_example.wrg", "w") as f:
+    # Write the top line of the file
+    f.write(f"{Nx} {Ny} {Xmin} {Ymin} {cell_size}\n")
+
+    # Now loop over the points
+    for i in range(Nx * Ny):
+        # Initiate the line to write as 10 blank spaces
+        line = "          "
+
+        # Add the x-coodinate as a 10 character fixed width integer
+        line = line + f"{int(x_list[i]):10d}"
+
+        # Add the y-coodinate as a 10 character fixed width integer
+        line = line + f"{int(y_list[i]):10d}"
+
+        # Add the z-coodinate as a 10 character fixed width integer
+        line = line + f"{int(z_coord):8d}"
+
+        # Add the height above ground level as a 10 character fixed width integer
+        line = line + f"{int(height_above_ground_level):5d}"
+
+        # Get the wind rose for this point
+        wind_rose = wind_rose_list[i]
+
+        # Get the frequency matrix and wind speed
+        freq_table = wind_rose.freq_table
+        wind_speeds = wind_rose.wind_speeds
+        wind_directions = wind_rose.wind_directions
+
+        # Get the A and k parameters across all sectors
+        freq_table_ws = freq_table.sum(axis=0)
+        A, k = estimate_weibull(wind_speeds, freq_table_ws)
+
+        # Write the A and k parameters
+        line = line + f"{A:5.1f}{k:6.2f}"
+
+        # Add place holder 0 for the power density
+        line = line + f"{0:15d}"
+
+        # Write the number of sectors
+        line = line + f"{num_sectors:3d}"
+
+        # Get the frequency table across wind directions
+        freq_table_wd = freq_table.sum(axis=1)
+
+        # Step through the sectors
+        for wd_idx in range(num_sectors):
+            # Write the probability for this sector
+            line = line + f"{int(1000*freq_table_wd[wd_idx]):4d}"
+
+            # Get the A and k parameters for this sector
+            A, k = estimate_weibull(wind_speeds, freq_table[wd_idx, :])
+
+            # Write the A and k parameters
+            line = line + f"{int(A*10):4d}{int(k*100):5d}"
+
+        # Write the line to the file
+        f.write(line + "\n")
+
+
+# Show the plots
+plt.show()
diff --git a/examples/examples_wind_resource_grid/001_wind_rose_wrg.py b/examples/examples_wind_resource_grid/001_wind_rose_wrg.py
@@ -0,0 +1,77 @@
+"""Example: WindRoseWRG
+
+`WindRoseWRG` is a type of WindData object, like `WindRose` and `TimeSeries`, that
+is used to store wind data in a format that can be used by the FLORIS model.  `WindRoseWRG`
+is different that `WindRose` however because the internal data holds the information
+of the WRG file and then a `WindRose` object is created for each turbine in a provided
+layout.
+
+In this example the WRG file generated in the previous example is read in
+using the `WindRoseWRG` object, and wind roses as points on the WRG grid, as will
+as in-between interpolated points have wind roses calculated using the `get_wind_rose_at_point`
+method.  Finally, the wind roses are upsampled to 5 degree wind direction bins and plotted.
+
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+from floris import WindRoseWRG
+
+
+# Read the WRG file
+wind_rose_wrg = WindRoseWRG("wrg_example.wrg")
+
+# Print some basic information
+print(wind_rose_wrg)
+
+# The wind roses were set to have a higher concentration of faster north winds for
+# increasing y, show that this is contained within the wind roses, even those interpolated
+# between grid points
+y_points_to_test = np.array([0, 500, 1000, 1500, 2000])
+
+fig, axarr = plt.subplots(1, 5, figsize=(16, 5), subplot_kw={"polar": True})
+
+for i in range(5):
+    wind_rose = wind_rose_wrg.get_wind_rose_at_point(0, y_points_to_test[i])
+    wind_rose.plot(ax=axarr[i], ws_step=5)
+    if i %2 == 0:
+        axarr[i].set_title(f"y = {y_points_to_test[i]}")
+    else:
+        axarr[i].set_title(f"y = {y_points_to_test[i]}\n(Interpolated)")
+
+# Go through the axarr and delete the legends except for the middle
+for ax in [axarr[0], axarr[1], axarr[3], axarr[4]]:
+    ax.legend().set_visible(False)
+
+
+# Draw a horizontal line on each axis indicating the level of the lower wind speed
+# bucket for the north wind from the first wind rose
+for i in range(5):
+    axarr[i].axhline(y=0.036, color="red", alpha=0.5)
+
+fig.suptitle("Wind Roses at locations with increasing y.  Note the location where the 5 m/s bin \
+transitions to 10 m/s for north wind at y = 0 is \nindicated by the red line to show \
+the increase in wind speed to the north as y increases.")
+
+# Since wind directions was not specified, the wind directions implied by the number of sectors
+# in the WRG was used, however the wind directions can be set using the set_wind_directions method
+# or passed in at initialization.  Here we upsample from 12, 30-deg sectors, to 72 5-deg sectors
+wind_rose_wrg.set_wd_step(5.0)
+
+fig, axarr = plt.subplots(1, 5, figsize=(16, 5), subplot_kw={"polar": True})
+
+for i in range(5):
+    wind_rose = wind_rose_wrg.get_wind_rose_at_point(0, y_points_to_test[i])
+    wind_rose.plot(ax=axarr[i], ws_step=5)
+    if i %2 == 0:
+        axarr[i].set_title(f"y = {y_points_to_test[i]}")
+    else:
+        axarr[i].set_title(f"y = {y_points_to_test[i]}\n(Interpolated)")
+
+# Go through the axarr and delete all the legends except for the middle
+for ax in axarr:
+    ax.legend().set_visible(False)
+
+fig.suptitle('Wind roses with upsampling to 5-deg bins')
+
+plt.show()