Add max_shading_elevation (#21)

* Squashed commit of the following:

commit f435e1c
Author: Adam R. Jensen <[email protected]>
Date:   Wed Mar 2 15:52:15 2022 +0100

    Add tests (#16)

    * Create twoaxistrackerfield.py

    * Add TwoAxisTrackerField to __init__.py and documentation.rst

    * Add _calculate_l_min to layout

    * Fix linter error

    * Update twoaxistrackerfield.py

    * Reorder layout ValueError warnings

    * Fix typo

    * Remove layout_type from generate_field_layout function

    * Add layout_type code to TwoAxisTrackerField

    * L_min -> min_tracker_spacing

    * Add whatsnew entry

    * Fix linting errors

    * Update twoaxistrackerfield.py

    * Correct horizon shading

    * Update twoaxistracking/twoaxistrackerfield.py

    Co-authored-by: Kevin Anderson <[email protected]>

    * Update twoaxistracking/twoaxistrackerfield.py

    Co-authored-by: Kevin Anderson <[email protected]>

    * Update twoaxistracking/twoaxistrackerfield.py

    Co-authored-by: Kevin Anderson <[email protected]>

    * Return scalar when scalar is passed

    * Correct use of np.isscalar

    * Update documentation

    * Update whatsnew

    * Rework intro_tutorial.ipynb

    * Update README.md

    * Create test_layout.py

    * Update test_layout.py

    * Use test folder in package folder

    * Delete test_placeholder.py

    * Add test files

    * Fix linting errors

    * Update scalar test

    * Add pandas to test in setup.cfg

    * Revert "package metadata and doc build updates (#17)"

    This reverts commit e9b5810.

    * Undo commits

    * Revert "Undo commits"

    This reverts commit 147ba2a.

    * Add test files

    * Switch to using np.testing.assert_allclose

    * Add test for sloped field layout

    * Move test_twoaxistracker to twoaxistracker PR

    * Add conftest.py to avoid duplicate fixtures

    * Remove imports of fixtures from conftest

    * Remove unnecessary imports

    * Increase test tolerance for test_field_slope

    * Add multi-polygon active area test

    * Add empty __init__.py

    * Use "from .conftest" as Kevin said..

    * Base square_field_layout_slope fixture off square_field_layout

    * Fix linter error

    * Remove backslash in conftest.py

    Co-authored-by: Kevin Anderson <[email protected]>

* Add get_max_shading_elevation func

Also, correct the relative slope calculation

* Correct relative_slope in conftest

* Add max elevation to TwoAxisTrackerField

* Update _field_layout_plot function

* Make get_horizon_elevation_angle separate function

* Add documentation entries

* Add test coverage

* Fix linter

* Minor documentation updates

* Update twoaxistracking_logo.svg

Include the shading from the base of the reference collector.

* Update twoaxistracking/layout.py

Co-authored-by: Kevin Anderson <[email protected]>

* Implement changes from review by kanderso-nrel

* Fix linter error

* Rewrite max_shading_elevation function

* Fix linter

* Add check for total area enclosing active areas

* Updates addressing review by kanderso-nrel

Co-authored-by: Kevin Anderson <[email protected]>

Author: AdamRJensen

docs/source/_static/twoaxistracking_logo.svg

@@ -12,3 +12,5 @@ Code documentation
    TwoAxisTrackerField
    TwoAxisTrackerField.get_shaded_fraction
    TwoAxisTrackerField.plot_field_layout
+   layout.max_shading_elevation
+   shading.horizon_elevation_angle

docs/source/documentation.rst

@@ -1,3 +1,4 @@
+# Import of functions that should be accessible from the package top-level
 from .layout import generate_field_layout  # noqa: F401
 from .shading import shaded_fraction  # noqa: F401
 from .twoaxistrackerfield import TwoAxisTrackerField  # noqa: F401

twoaxistracking/__init__.py

@@ -121,6 +121,65 @@ def generate_field_layout(gcr, total_collector_area, min_tracker_spacing,
     relative_azimuth = np.mod(450-np.rad2deg(np.arctan2(Y, X)), 360)
     # Relative slope of collectors
     # positive means collector is higher than reference collector
-    relative_slope = -np.cos(np.deg2rad(slope_azimuth - relative_azimuth)) * slope_tilt  # noqa: E501
+    relative_slope = np.rad2deg(np.arctan(-np.cos(np.deg2rad(slope_azimuth - relative_azimuth))
+                                          * np.tan(np.deg2rad(slope_tilt))))
 
     return X, Y, Z, tracker_distance, relative_azimuth, relative_slope
+
+
+def max_shading_elevation(total_collector_geometry, tracker_distance,
+                          relative_slope):
+    """Calculate the maximum elevation angle for which shading can occur.
+
+    Parameters
+    ----------
+    total_collector_geometry: Shapely Polygon
+        Polygon corresponding to the total collector area.
+    tracker_distance: array of floats
+        Distances between neighboring trackers and the reference tracker.
+    relative_slope: array of floats
+        Slope between neighboring trackers and reference tracker. A positive
+        slope means neighboring collector is higher than reference collector.
+
+    Returns
+    -------
+    max_shading_elevation: float
+        The highest solar elevation angle for which shading can occur for a
+        given field layout and collector geometry [degrees]
+
+    Note
+    ----
+    The maximum shading elevation angle is calculated for all neighboring
+    trackers using the bounding box geometry and the bounding circle. For
+    rectangular collectors (as approximated when using the bounding box), the
+    maximum shading elevation occurs when one of the upper corners of the
+    projected shading geometry and the lower corner of the reference collector
+    intersects. For circular collectors (as approximated by the bounding
+    cirlce), the maximum elevation occurs when the projected shadow is directly
+    below the reference collector and the two circles tangent to each other.
+
+    The maximum elevation is calculated using both the bounding box and the
+    bounding circle, and the minimum of these two elevations is returned. For
+    rectangular and circular collectors, the maximum elevation is exact,
+    whereas for other geometries, the returned elevation is a conservative
+    estimate.
+    """
+    # Calculate extent of box bounding the total collector geometry
+    x_min, y_min, x_max, y_max = total_collector_geometry.bounds
+    # Collector dimensions
+    x_dim = x_max - x_min
+    y_dim = y_max - y_min
+    delta_gamma_rad = np.arcsin(x_dim / tracker_distance)
+    # Calculate max elevation based on the bounding box (rectangular)
+    max_elevations_rectangular = np.rad2deg(np.arcsin(
+        y_dim * np.cos(np.deg2rad(relative_slope)) /
+        (tracker_distance * np.cos(delta_gamma_rad)))) + relative_slope
+    # Calculate max elevations using the minimum bounding diameter (circular)
+    D_min = _calculate_min_tracker_spacing(total_collector_geometry)
+    max_elevations_circular = np.rad2deg(np.arcsin(
+        (D_min * np.cos(np.deg2rad(relative_slope)))/tracker_distance)) \
+        + relative_slope
+    # Compute max elevation
+    max_elevation = np.nanmin([np.nanmax(max_elevations_rectangular),
+                               np.nanmax(max_elevations_circular)])
+    return max_elevation

twoaxistracking/layout.py

@@ -15,7 +15,7 @@ def _plot_field_layout(X, Y, Z, min_tracker_spacing):
     # to correctly display the middle color when all tracker Z coords are zero
     cmap = cm.viridis_r
     colors = cmap(norm(Z))
-    fig, ax = plt.subplots(figsize=(6, 6), subplot_kw={'aspect': 'equal'})
+    fig, ax = plt.subplots(figsize=(4, 4), subplot_kw={'aspect': 'equal'})
     # Plot a circle for each neighboring collector (diameter equals min_tracker_spacing)
     ax.add_collection(collections.EllipseCollection(
         widths=min_tracker_spacing, heights=min_tracker_spacing, angles=0,
@@ -26,6 +26,8 @@ def _plot_field_layout(X, Y, Z, min_tracker_spacing):
         widths=min_tracker_spacing, heights=min_tracker_spacing, angles=0,
         units='xy', facecolors='red', edgecolors=("black",), linewidths=(1,),
         offsets=[0, 0], transOffset=ax.transData))
+    ax.set_xlabel('Tracker position (east-west direction)')
+    ax.set_ylabel('Tracker position (north-south direction)')
     fig.colorbar(cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax, shrink=0.8,
                  label='Relative tracker height (vertical)')
     # Set limits
@@ -62,7 +64,7 @@ def _plot_shading(active_collector_geometry, unshaded_geometry,
         shading_geometries, facecolor='blue', linewidth=0.5, alpha=0.5)
 
     fig, axes = plt.subplots(1, 2, subplot_kw=dict(aspect='equal'))
-    axes[0].set_title('Total area and shading areas')
+    axes[0].set_title('Active area and shading areas')
     axes[0].add_collection(active_patches, autolim=True)
     axes[0].add_collection(shading_patches, autolim=True)
     axes[1].set_title('Unshaded area')

twoaxistracking/plotting.py

@@ -3,10 +3,37 @@
 from twoaxistracking import plotting
 
 
+def horizon_elevation_angle(azimuth, slope_azimuth, slope_tilt):
+    """Calculate horizon elevation angle caused by a sloped field.
+
+    Parameters
+    ----------
+    azimuth : float
+        Azimuth angle for which the horizon elevation angle is to be calculated
+        [degrees]
+    slope_azimuth : float
+        Direction of normal to slope on horizontal [degrees]
+    slope_tilt : float
+        Tilt of slope relative to horizontal [degrees]
+
+    Returns
+    -------
+    horizon_elevation_angle : float
+        Horizon elevation angle [degrees]
+    """
+    horizon_elevation_angle = np.rad2deg(np.arctan(
+        - np.cos(np.deg2rad(slope_azimuth - azimuth))
+        * np.tan(np.deg2rad(slope_tilt))))
+    # Horizon elevation angle cannot be less than zero
+    horizon_elevation_angle = np.clip(horizon_elevation_angle, a_min=0, a_max=None)
+    return horizon_elevation_angle
+
+
 def shaded_fraction(solar_elevation, solar_azimuth,
                     total_collector_geometry, active_collector_geometry,
                     min_tracker_spacing, tracker_distance, relative_azimuth,
-                    relative_slope, slope_azimuth=0, slope_tilt=0, plot=False):
+                    relative_slope, slope_azimuth=0, slope_tilt=0,
+                    max_shading_elevation=90, plot=False):
     """Calculate the shaded fraction for any layout of two-axis tracking collectors.
 
     Parameters
@@ -29,12 +56,18 @@ def shaded_fraction(solar_elevation, solar_azimuth,
     relative_slope: array of floats
         Slope between neighboring trackers and reference tracker. A positive
         slope means neighboring collector is higher than reference collector.
-    slope_azimuth : float
+    slope_azimuth : float, optional
         Direction of normal to slope on horizontal [degrees]. Used to determine
         horizon shading.
-    slope_tilt : float
+    slope_tilt : float, optional
         Tilt of slope relative to horizontal [degrees]. Used to determine
         horizon shading.
+    max_shading_elevation : float, optional
+        The maximum elevation angle for which shading may occur. Specifying the
+        max_shading_elevation skips the calculations for which the solar
+        elevation angle is higher and sets the shaded fraction to zero.
+        This reduces the calculation time and results in the same output
+        assuming the correct value has been provided.
     plot: bool, default: True
         Whether to plot the projected shadows and unshaded area.
 
@@ -46,11 +79,13 @@ def shaded_fraction(solar_elevation, solar_azimuth,
     # If the sun is below the horizon, set the shaded fraction to nan
     if solar_elevation < 0:
         return np.nan
+    # Set shaded fraction to 0 (unshaded) if solar elevation is higher than
+    # max_shading_elevation
+    elif solar_elevation > max_shading_elevation:
+        return 0
     # Set shaded fraction to 1 (fully shaded) if the solar elevation is below
     # the horizon line caused by the tilted ground
-    elif np.tan(np.deg2rad(solar_elevation)) <= (
-            - np.cos(np.deg2rad(slope_azimuth-solar_azimuth))
-            * np.tan(np.deg2rad(slope_tilt))):
+    elif solar_elevation <= horizon_elevation_angle(solar_azimuth, slope_azimuth, slope_tilt):
         return 1
 
     azimuth_difference = solar_azimuth - relative_azimuth

twoaxistracking/shading.py

@@ -12,6 +12,15 @@ def rectangular_geometry():
     return collector_geometry, total_collector_area, min_tracker_spacing
 
 
+@pytest.fixture
+def circular_geometry():
+    # A circular collector centered at (0,0) and has a radius of 2
+    collector_geometry = geometry.Point(0, 0).buffer(2)
+    total_collector_area = collector_geometry.area
+    min_tracker_spacing = layout._calculate_min_tracker_spacing(collector_geometry)
+    return collector_geometry, total_collector_area, min_tracker_spacing
+
+
 @pytest.fixture
 def active_geometry_split():
     active_collector_geometry = geometry.MultiPolygon([
@@ -41,8 +50,8 @@ def square_field_layout_sloped(square_field_layout):
     X, Y, _, tracker_distance, relative_azimuth, _ = square_field_layout
     Z = np.array([0.12372765, 0.06186383, 0, 0.06186383,
                   -0.06186383, 0, -0.06186383, -0.12372765])
-    relative_slope = np.array([5, 3.53553391, 0, 3.53553391,
-                               -3.53553391, 0, -3.53553391, -5])
+    relative_slope = np.array([5, 3.540025, 0, 3.540025,
+                               -3.540025, 0, -3.540025, -5])
     return X, Y, Z, tracker_distance, relative_azimuth, relative_slope
 
 

twoaxistracking/tests/conftest.py

@@ -145,3 +145,33 @@ def test_neighbor_order(rectangular_geometry):
             offset=0,
             rotation=0)
     assert len(X) == (7*7-1)
+
+
+def test_calculation_of_max_shading_elevation_rectangle(rectangular_geometry, square_field_layout):
+    # Test that the maximum elevation angle for which shading can occur is
+    # calculated correctly for rectangular collectors
+    collector_geometry, total_collector_area, min_tracker_spacing = rectangular_geometry
+    X, Y, Z, tracker_distance, relative_azimuth, relative_slope = \
+        square_field_layout
+    max_shading_elevation = layout.max_shading_elevation(
+        collector_geometry, tracker_distance, relative_slope)
+    np.testing.assert_allclose(max_shading_elevation, 16.77865488)
+
+
+def test_calculation_of_max_shading_elevation_circle(circular_geometry):
+    # Test that the maximum elevation angle for which shading can occur is
+    # calculated correctly for closely packed circular collectors
+    collector_geometry, total_collector_area, min_tracker_spacing = circular_geometry
+    X, Y, Z, tracker_distance, relative_azimuth, relative_slope = \
+        layout.generate_field_layout(
+            gcr=0.5,
+            total_collector_area=total_collector_area,
+            min_tracker_spacing=min_tracker_spacing,
+            neighbor_order=2,
+            aspect_ratio=1,
+            offset=0,
+            rotation=0)
+
+    max_shading_elevation = layout.max_shading_elevation(
+        collector_geometry, tracker_distance, relative_slope)
+    np.testing.assert_allclose(max_shading_elevation, 43.16784217)

twoaxistracking/tests/test_layout.py

@@ -101,3 +101,25 @@ def test_shading_below_hill_horizon(rectangular_geometry, square_field_layout):
         slope_tilt=10,
         plot=False)
     assert shaded_fraction == 1
+
+
+def test_shading_max_shading_elevation(rectangular_geometry, square_field_layout):
+    # Test that shaded_fraction is set to one when the solar elevation angle
+    # is greater than the max_shading_elevation (even though shading may occur)
+    collector_geometry, total_collector_area, min_tracker_spacing = rectangular_geometry
+    X, Y, Z, tracker_distance, relative_azimuth, relative_slope = \
+        square_field_layout
+    shaded_fraction = shading.shaded_fraction(
+        solar_elevation=3,  # low solar elevation angle with guaranteed shading
+        solar_azimuth=180,
+        total_collector_geometry=collector_geometry,
+        active_collector_geometry=collector_geometry,
+        min_tracker_spacing=min_tracker_spacing,
+        tracker_distance=tracker_distance,
+        relative_azimuth=relative_azimuth,
+        relative_slope=relative_slope,
+        slope_azimuth=0,
+        slope_tilt=10,
+        max_shading_elevation=2,  # lower than true max angle for testing purposes
+        plot=False)
+    assert shaded_fraction == 0

twoaxistracking/tests/test_shading.py

@@ -179,6 +179,8 @@ def test_calculation_of_shaded_fraction_array(rectangular_geometry, solar_positi
 def test_calculation_of_shaded_fraction_float(rectangular_geometry):
     # Test if shaded fraction is calculated correct when solar elevation and
     # azimuth are scalar
+    # Also tests that no error is raised when total and active geometries are
+    # identical.
     collector_geometry, total_collector_area, min_tracker_spacing = rectangular_geometry
     field = twoaxistrackerfield.TwoAxisTrackerField(
         total_collector_geometry=collector_geometry,
@@ -192,3 +194,19 @@ def test_calculation_of_shaded_fraction_float(rectangular_geometry):
     result = field.get_shaded_fraction(40, 180)
     np.testing.assert_allclose(result, 0)
     assert np.isscalar(result)
+
+
+def test_total_collector_geometry_encloses_active_areas(rectangular_geometry, circular_geometry):
+    # Test that ValueError is raised if the aperture collector geometry is not
+    # completely enclosed by the total collector geometry
+    rectangular_collector, total_collector_area, min_tracker_spacing = rectangular_geometry
+    circular_collector, total_collector_area, min_tracker_spacing = circular_geometry
+    with pytest.raises(ValueError, match="does not completely enclose"):
+        _ = twoaxistrackerfield.TwoAxisTrackerField(
+            total_collector_geometry=rectangular_collector,
+            active_collector_geometry=circular_collector,
+            neighbor_order=1,
+            gcr=0.1,
+            aspect_ratio=1,
+            offset=0,
+            rotation=0)