diff --git a/.gitignore b/.gitignore index e124a08fdb..1605e1780f 100644 --- a/.gitignore +++ b/.gitignore @@ -97,4 +97,3 @@ coverage.xml # airspeed velocity env results - diff --git a/ci/requirements-py3.10.yml b/ci/requirements-py3.10.yml index fcb89d2e7f..89cb593ff7 100644 --- a/ci/requirements-py3.10.yml +++ b/ci/requirements-py3.10.yml @@ -22,7 +22,9 @@ dependencies: - pytz - requests - scipy >= 1.6.0 + - shapely >= 2.0.2 - statsmodels - pip: - nrel-pysam>=2.0 - solarfactors + - shapely >= 2.0.2 diff --git a/ci/requirements-py3.11.yml b/ci/requirements-py3.11.yml index f6556ecf94..5be297e6ab 100644 --- a/ci/requirements-py3.11.yml +++ b/ci/requirements-py3.11.yml @@ -22,7 +22,9 @@ dependencies: - pytz - requests - scipy >= 1.6.0 + - shapely >= 2.0.2 - statsmodels - pip: - nrel-pysam>=2.0 - solarfactors + - shapely >= 2.0.2 diff --git a/ci/requirements-py3.12.yml b/ci/requirements-py3.12.yml index f3d8fc2d0c..489350aac0 100644 --- a/ci/requirements-py3.12.yml +++ b/ci/requirements-py3.12.yml @@ -22,7 +22,9 @@ dependencies: - pytz - requests - scipy >= 1.6.0 + - shapely >= 2.0.2 - statsmodels - pip: - nrel-pysam>=2.0 # - solarfactors # required shapely<2 isn't available for 3.12 + - shapely >= 2.0.2 diff --git a/ci/requirements-py3.8-min.yml b/ci/requirements-py3.8-min.yml index ef230b1627..aed65ddb32 100644 --- a/ci/requirements-py3.8-min.yml +++ b/ci/requirements-py3.8-min.yml @@ -19,4 +19,5 @@ dependencies: - scipy==1.6.0 - pytest-rerunfailures # conda version is >3.6 - pytest-remotedata # conda package is 0.3.0, needs > 0.3.1 + - shapely >= 2.0.2 - requests-mock diff --git a/ci/requirements-py3.8.yml b/ci/requirements-py3.8.yml index 19c3b8f2f3..0ed6bfa119 100644 --- a/ci/requirements-py3.8.yml +++ b/ci/requirements-py3.8.yml @@ -26,3 +26,4 @@ dependencies: - pip: - nrel-pysam>=2.0 - solarfactors + - shapely >= 2.0.2 diff --git a/ci/requirements-py3.9.yml b/ci/requirements-py3.9.yml index b5aa976b4b..19a7194b7d 100644 --- a/ci/requirements-py3.9.yml +++ b/ci/requirements-py3.9.yml @@ -22,7 +22,9 @@ dependencies: - pytz - requests - scipy >= 1.6.0 + - shapely >= 2.0.2 - statsmodels - pip: - nrel-pysam>=2.0 - - solarfactors \ No newline at end of file + - solarfactors + - shapely >= 2.0.2 diff --git a/docs/examples/agrivoltaics/README.rst b/docs/examples/agrivoltaics/README.rst new file mode 100644 index 0000000000..9c3c2ad045 --- /dev/null +++ b/docs/examples/agrivoltaics/README.rst @@ -0,0 +1,2 @@ +Agrivoltaic Systems Modelling +----------------------------- diff --git a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py new file mode 100644 index 0000000000..c69b1fc8b4 --- /dev/null +++ b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py @@ -0,0 +1,291 @@ +""" +Fixed Tilt System impact on crop shading +======================================== + +This example shows how to calculate the shading of a crop field by a fixed tilt +system. +""" + +# %% +# This is the first of a series of examples that will show how to calculate the +# shading of a crop field by a fixed tilt system, a single-axis tracker, and a +# two-axis tracker. The examples will show how to calculate the shading in 3D +# and 2D space, and how to calculate the shading fraction with the help of the +# :py:mod:`~pvlib.spatial` module and its classes. +# +# Paper Examples +# ============== +# +# +--------------------+----------+-------------+----------+-------+ +# | Input parameters | Vertical | Single-axis | Two-axis | Units | +# +====================+==========+=============+==========+=======+ +# | Panel width | 1 | 1 | 1 | [m] | +# +--------------------+----------+-------------+----------+-------+ +# | Panel length | 2 | 2 | 2 | [m] | +# +--------------------+----------+-------------+----------+-------+ +# | Number of panels | 40 | 40 | 40 | [-] | +# +--------------------+----------+-------------+----------+-------+ +# | Total panel area | 80 | 80 | 80 | [m²] | +# +--------------------+----------+-------------+----------+-------+ +# | Number of rows | 2 | 2 | 2 | [-] | +# +--------------------+----------+-------------+----------+-------+ +# | Row spacing | 10 | 10 | 10 | [m] | +# +--------------------+----------+-------------+----------+-------+ +# | Row length | 20 | 20 | 20 | [m] | +# +--------------------+----------+-------------+----------+-------+ +# | Crop area | 200 | 200 | 200 | [m²] | +# +--------------------+----------+-------------+----------+-------+ +# | Pitch | \- | \- | 2 | [m] | +# +--------------------+----------+-------------+----------+-------+ +# | Height | 0 | 3 | 3 | [m] | +# +--------------------+----------+-------------+----------+-------+ +# | Fixed tilt angle | 90 | \- | \- | [°] | +# +--------------------+----------+-------------+----------+-------+ +# | Azimuth angle | 0 | 0 | 0 | [°] | +# +--------------------+----------+-------------+----------+-------+ +# | Maximum tilt angle | \- | 60 | 60 | [°] | +# +--------------------+----------+-------------+----------+-------+ +# | Minimum tilt angle | \- | -60 | -60 | [°] | +# +--------------------+----------+-------------+----------+-------+ +# +# + +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d.art3d import Poly3DCollection +import matplotlib.animation as animation +from matplotlib.dates import DateFormatter +import shapely +try: + from shapely import Polygon # shapely >= 2 +except ImportError: + from shapely.geometry import Polygon # shapely < 2 +import pandas as pd +import numpy as np +from functools import partial +from pvlib.spatial import RectangularSurface +from pvlib import solarposition + +# Kärrbo Prästgård, Västerås, Sweden +latitude, longitude, altitude = 59.6099, 16.5448, 20 # °N, °E, m + +spring_equinox = pd.date_range("2021-03-20", periods=24, freq="H") +summer_solstice = pd.date_range("2021-06-21", periods=24, freq="H") +fall_equinox = pd.date_range("2021-09-22", periods=24, freq="H") +winter_solstice = pd.date_range("2021-12-21", periods=24, freq="H") +dates = ( + spring_equinox.union(summer_solstice) + .union(fall_equinox) + .union(winter_solstice) +) +# dates = spring_equinox +solar_position = solarposition.get_solarposition( + dates, latitude, longitude, altitude +) +solar_zenith = solar_position["apparent_zenith"] +solar_azimuth = solar_position["azimuth"] +N = len(solar_zenith) + +# %% +# Fixed Tilt +# ---------- +# The fixed tilt system is composed of a crop field and two vertical rows of +# PV panels. Rows are placed at the long sides of the field, with the panels +# facing east and west. The field is 20 m long and 10 m wide, and the panels +# are 2 m wide (height since they are vertical) and 20 m long. + +field = RectangularSurface( # crops surface + center=[5, 10, 0], + azimuth=90, + tilt=0, + axis_tilt=0, + width=10, + length=20, +) +pv_rows = ( + RectangularSurface( # west-most row (lowest X-coordinate) + center=[-10 / 2 + 5, 10, 2 / 2], + azimuth=90, + tilt=90, + axis_tilt=0, + width=2, + length=20, + ), + RectangularSurface( # east-most row (highest X-coordinate) + center=[10 / 2 + 5, 10, 2 / 2], + azimuth=90, + tilt=90, + axis_tilt=0, + width=2, + length=20, + ), +) + + +# %% +# Run the simulation +# ------------------ +# The shading fraction is calculated at each instant in time, and the results +# are stored in the `fixed_tilt_shaded_fraction` array. This is done because +# the shading calculation API does not allow for vectorized calculations. + +# Allocate space for the shading results +fixed_tilt_shaded_fraction = np.zeros((N,), dtype=float) + + +# Shades callback +def simulation_and_plot_callback( + timestamp_index, *, shade_3d_artists, shade_2d_artists, sun_annotation +): + # Calculate the shades at an specific instant in time + solar_zenith_instant = solar_zenith.iloc[timestamp_index] + solar_azimuth_instant = solar_azimuth.iloc[timestamp_index] + # Update the sun position text + sun_annotation.set_text( + f"Sun at γ={solar_zenith_instant:.2f}, β={solar_azimuth_instant:.2f}" + ) + # skip this instant if the sun is below the horizon + if solar_zenith_instant < 0: + fixed_tilt_shaded_fraction[timestamp_index] = 0 + return *shade_3d_artists, *shade_2d_artists + # Calculate the shades, both in 3D and 2D + shades_3d = field.get_3D_shades_from( + solar_zenith_instant, solar_azimuth_instant, *pv_rows + ) + shades_2d = field.get_2D_shades_from( + solar_zenith_instant, solar_azimuth_instant, shades_3d=shades_3d + ) + # Plot the calculated shades + for index, shade in enumerate(shades_3d.geoms): # 3D + if not shade.is_empty: + shade_3d_artists[index].set_verts( + [shade.exterior.coords], + closed=False, # polygon is already closed + ) + for index, shade in enumerate(shades_2d.geoms): # 2D + if not shade.is_empty: + shade_2d_artists[index].set_path( + shapely.plotting._path_from_polygon(shade) + ) + + # Calculate the shaded fraction + fixed_tilt_shaded_fraction[timestamp_index] = ( + sum(shade.area for shade in shades_2d.geoms) / field2d.area + ) + return *shade_3d_artists, *shade_2d_artists + + +# %% +# Plot both the 3D and 2D shades +# ------------------------------ + +fig = plt.figure(figsize=(10, 10)) +gs = fig.add_gridspec(10, 1) + +# Plotting styles +field_style = {"color": "forestgreen", "alpha": 0.5} +row_style = {"color": "darkblue", "alpha": 0.5} +shade_style = {"color": "dimgrey", "alpha": 0.8} +field_style_2d = {**field_style, "add_points": False} +shade_style_2d = {**shade_style, "add_points": False} + +ax1 = fig.add_subplot(gs[0:6, 0], projection="3d") +ax2 = fig.add_subplot(gs[8:, 0]) + +# Upper plot, 3D +ax1.axis("equal") +ax1.view_init( + elev=45, + azim=-60, # matplotlib's azimuth is right-handed to Z+, measured from X+ +) +ax1.set_xlim(-0, 15) +ax1.set_ylim(0, 20) +ax1.set_zlim(0, 10) +ax1.set_xlabel("West(-) East(+) [m]") +ax1.set_ylabel("South(-) North(+) [m]") +field.plot(ax=ax1, **field_style) +for pv_row in pv_rows: + pv_row.plot(ax=ax1, **row_style) + +# Lower plot, 2D +field2d = field.representation_in_2D_space() +shapely.plotting.plot_polygon(field2d, ax=ax2, **field_style_2d) + +# Add empty shade artists for each shading object, in this case each of the +# PV rows. Artists will be updated in the animation callback later. +shade3d_artists = ( + ax1.add_collection3d(Poly3DCollection([], **shade_style)), +) * len(pv_rows) +shade2d_artists = ( + shapely.plotting.plot_polygon( + Polygon([[0, 0]] * 4), ax=ax2, **shade_style_2d + ), +) * len(pv_rows) +sun_text_artist = fig.text(0.5, 0.95, "Sun at γ=--, β=--", ha="center") + +ani = animation.FuncAnimation( + fig, + partial( + simulation_and_plot_callback, + shade_3d_artists=shade3d_artists, + shade_2d_artists=shade2d_artists, + sun_annotation=sun_text_artist, + ), + frames=np.arange(N), + interval=200, + blit=True, +) + +# uncomment to run and save animation locally +# ani.save("fixed_tilt_shading.gif", writer="pillow") + +plt.show() + +# %% +# Shaded Fraction vs. Time +# ------------------------ +# Create a handy pandas series to plot the shaded fraction vs. time. + +fixed_tilt_shaded_fraction = pd.Series(fixed_tilt_shaded_fraction, index=dates) + +fig, axs = plt.subplots(ncols=4, sharey=True, figsize=(20, 5)) +fig.suptitle("Shaded Fraction vs. Time") +fig.subplots_adjust(wspace=0) + +for ax, day_datetimes, title in zip( + axs, + (spring_equinox, summer_solstice, fall_equinox, winter_solstice), + ("Spring Equinox", "Summer Solstice", "Autumn Equinox", "Winter Solstice"), +): + fixed_tilt_shaded_fraction[day_datetimes].plot(ax=ax) + ax.xaxis.set_major_formatter(DateFormatter("%H")) + ax.set_title(title) + ax.grid(True) +for ax_a, ax_b in zip(axs[:-1], axs[1:]): + ax_a.spines.right.set_visible(False) + ax_b.spines.left.set_visible(False) + ax_a.tick_params(labelright=False) + ax_b.tick_params(labelleft=False) + +axs[0].set_ylabel("Shaded Fraction [Unitless]") +axs[0].set_ylim(0, 1) + +plt.show() + + +# %% +# References +# ---------- +# .. [1] S. Zainali et al., 'Direct and diffuse shading factors modelling +# for the most representative agrivoltaic system layouts', Applied +# Energy, vol. 339, p. 120981, Jun. 2023, +# :doi:`10.1016/j.apenergy.2023.120981`. +# .. [2] Y. Cascone, V. Corrado, and V. Serra, 'Calculation procedure of +# the shading factor under complex boundary conditions', Solar Energy, +# vol. 85, no. 10, pp. 2524-2539, Oct. 2011, +# :doi:`10.1016/j.solener.2011.07.011`. +# .. [3] Kevin S. Anderson, Adam R. Jensen; Shaded fraction and +# backtracking +# in single-axis trackers on rolling terrain. J. Renewable Sustainable +# Energy 1 March 2024; 16 (2): 023504. :doi:`10.1063/5.0202220`. + +# %% diff --git a/docs/examples/shading/plot_spatial_row_to_row_shading.py b/docs/examples/shading/plot_spatial_row_to_row_shading.py new file mode 100644 index 0000000000..9fa77d8079 --- /dev/null +++ b/docs/examples/shading/plot_spatial_row_to_row_shading.py @@ -0,0 +1,115 @@ +""" +Spatial shading of a row-to-row shading system +============================================== +""" + +# %% +# This example demonstrates how to calculate the shading between two rows of +# panels in a row-to-row shading system. The example will show how to calculate +# the shaded fraction in 3D and 2D space with the help of the +# :py:mod:`~pvlib.spatial` module and its classes. +# +# This is a basic example on how to calculate and plot the shaded fraction +# for an instantaneous time. A more complex task is to calculate the shadows +# for a time range. This can be done by iterating over the time range and +# calculating the shadows for each time step. This is done since the +# :py:class:`~pvlib.spatial.FlatSurface` does not support the calculation of +# the shaded fraction for a time range. +# The example :ref:`sphx_glr_gallery_plot_par_direct_shading0_fixed_tilt.py` +# shows how to calculate the shading for a time range for a fixed tilt system. + +from pvlib.spatial import RectangularSurface +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d.art3d import Poly3DCollection +import shapely + +# %% +# Description of the system +# ------------------------- +# Let's start by creating a row-to-row system. We will create a +# rectangular surface for each of the two rows of panels. Both rows will be +# parallel to each other and will have the same azimuth angle. The tilt angle +# of the first row will be 20 degrees, and the tilt angle of the second row +# will be 30 degrees. The distance between the two rows will be 3 meters. +# Also, we will assume scalar values for the solar azimuth and zenith angles. +# Feel free to download the example and change the values to see how the +# shades change. + +solar_azimuth = 165 # degrees +solar_zenith = 75 # degrees + +row1 = RectangularSurface( # south-most row + center=[0, 0, 3], azimuth=165, tilt=20, axis_tilt=10, width=2, length=20 +) + +row2 = RectangularSurface( # north-most row + center=[0, 3, 3], azimuth=165, tilt=30, axis_tilt=10, width=2, length=20 +) + +# %% +# Calculating the shadows +# ----------------------- +# The 3D shapely polygons representing the shadows can be calculated with the +# :py:meth:`~pvlib.spatial.RectangularSurface.get_3D_shades_from` method. +# The 2D shapely polygons representing the shadows can be calculated with the +# :py:meth:`~pvlib.spatial.RectangularSurface.get_2D_shades_from` method. This +# method uses either the 3D shadows or calculates them internally if not +# provided. If the 3D shadows are needed outside its scope, it is recommended +# to calculate them separately and pass them as an argument for performance +# reasons. + +shades_3d = row2.get_3D_shades_from(solar_zenith, solar_azimuth, row1) +shades_2d = row2.get_2D_shades_from( + solar_zenith, solar_azimuth, shades_3d=shades_3d +) + +# %% +# Scene and shades plot +# --------------------- + +row_style = {"color": "darkblue", "alpha": 0.5} +shade_style = {"color": "dimgrey", "alpha": 0.8} +row_style_2d = {**row_style, "add_points": False} +shade_style_2d = {**shade_style, "add_points": False} + +fig = plt.figure(figsize=(10, 10)) + +# Split the figure in two axes +gs = fig.add_gridspec(10, 1) +ax1 = fig.add_subplot(gs[0:7, 0], projection="3d") +ax2 = fig.add_subplot(gs[8:, 0]) + +# 3D plot +ax1.view_init( + elev=60, + azim=-30, # matplotlib's azimuth is right-handed to Z+, measured from X+ +) +row1.plot(ax=ax1, **row_style) +row2.plot(ax=ax1, **row_style) +for shade in shades_3d.geoms: + if shade.is_empty: + continue # skip empty shades; else an exception will be raised + # use Matplotlib's Poly3DCollection natively since experimental + # shapely.plotting.plot_polygon does not support 3D + vertexes = shade.exterior.coords[:-1] + ax1.add_collection3d(Poly3DCollection([vertexes], **shade_style)) + +ax1.axis("equal") +ax1.set_zlim(0) +ax1.set_xlabel("West(-) East(+) [m]") +ax1.set_ylabel("South(-) North(+) [m]") + +# 2D plot +row2_2d = row2.representation_in_2D_space() +shapely.plotting.plot_polygon(row2_2d, ax=ax2, **row_style_2d) +for shade in shades_2d.geoms: + shapely.plotting.plot_polygon(shade, ax=ax2, **shade_style_2d) + +# %% +# Calculate the shaded fraction +# ----------------------------- +# The shaded fraction can be calculated by dividing the sum of the areas of the +# shadows by the area of the surface. The shaded fraction is a scalar value. + +shaded_fraction = sum(shade.area for shade in shades_2d.geoms) / row2_2d.area +print(f"The shaded fraction is {shaded_fraction:.2f}") diff --git a/docs/sphinx/source/conf.py b/docs/sphinx/source/conf.py index 9c0d90990e..9f44c0518f 100644 --- a/docs/sphinx/source/conf.py +++ b/docs/sphinx/source/conf.py @@ -372,6 +372,9 @@ def setup(app): # https://sphinx-gallery.github.io/dev/configuration.html#removing-config-comments # noqa: E501 'remove_config_comments': True, + + # allow creating animations with Matplotlib + 'matplotlib_animations': True, } # supress warnings in gallery output # https://sphinx-gallery.github.io/stable/configuration.html diff --git a/docs/sphinx/source/reference/index.rst b/docs/sphinx/source/reference/index.rst index 9083f85bdd..00809c65a7 100644 --- a/docs/sphinx/source/reference/index.rst +++ b/docs/sphinx/source/reference/index.rst @@ -20,3 +20,4 @@ API reference bifacial scaling location + spatial diff --git a/docs/sphinx/source/reference/spatial.rst b/docs/sphinx/source/reference/spatial.rst new file mode 100644 index 0000000000..13078401c0 --- /dev/null +++ b/docs/sphinx/source/reference/spatial.rst @@ -0,0 +1,12 @@ +.. currentmodule:: pvlib + +Spatial +======= + +Objects for creating 3D scenes and calculating shades + +.. autosummary:: + :toctree: generated/ + + spatial.FlatSurface + spatial.RectangularSurface diff --git a/pvlib/__init__.py b/pvlib/__init__.py index b5b07866a4..58e064b39f 100644 --- a/pvlib/__init__.py +++ b/pvlib/__init__.py @@ -24,6 +24,7 @@ soiling, solarposition, spa, + spatial, temperature, tools, tracking, diff --git a/pvlib/spatial.py b/pvlib/spatial.py new file mode 100644 index 0000000000..08f6e9e95c --- /dev/null +++ b/pvlib/spatial.py @@ -0,0 +1,720 @@ +""" +Spatial functions for shading and 3D scenes analysis. +""" + +from pvlib.tools import sind, cosd +import numpy as np +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + +import shapely + +try: + # import objects from shapely >= 2 + from shapely import Polygon, MultiPolygon +except ImportError: + # import objects from shapely < 2 + from shapely.geometry import Polygon, MultiPolygon + +if not hasattr(shapely, "Polygon"): + shapely.Polygon = Polygon + shapely.MultiPolygon = MultiPolygon + + +# --- BEGIN backport of shapely plotting capabilities --- +def _default_ax(): + import matplotlib.pyplot as plt + + ax = plt.gca() + ax.grid(True) + ax.set_aspect("equal") + return ax + + +def _path_from_polygon(polygon): + from matplotlib.path import Path + + if isinstance(polygon, shapely.MultiPolygon): + return Path.make_compound_path( + *[_path_from_polygon(poly) for poly in polygon.geoms] + ) + else: + return Path.make_compound_path( + Path(np.asarray(polygon.exterior.coords)[:, :2]), + *[ + Path(np.asarray(ring.coords)[:, :2]) + for ring in polygon.interiors + ], + ) + + +def patch_from_polygon(polygon, **kwargs): + """ + Gets a Matplotlib patch from a (Multi)Polygon. + + Note: this function is experimental, and mainly targeting (interactive) + exploration, debugging and illustration purposes. + + Parameters + ---------- + polygon : shapely.Polygon or shapely.MultiPolygon + **kwargs + Additional keyword arguments passed to the matplotlib Patch. + + Returns + ------- + Matplotlib artist (PathPatch) + """ + from matplotlib.patches import PathPatch + + return PathPatch(_path_from_polygon(polygon), **kwargs) + + +def plot_polygon( + polygon, + ax=None, + add_points=True, + color=None, + facecolor=None, + edgecolor=None, + linewidth=None, + **kwargs, +): + """ + Plot a (Multi)Polygon. + + Note: this function is experimental, and mainly targeting (interactive) + exploration, debugging and illustration purposes. + + Parameters + ---------- + polygon : shapely.Polygon or shapely.MultiPolygon + ax : matplotlib Axes, default None + The axes on which to draw the plot. If not specified, will get the + current active axes or create a new figure. + add_points : bool, default True + If True, also plot the coordinates (vertices) as points. + color : matplotlib color specification + Color for both the polygon fill (face) and boundary (edge). By default, + the fill is using an alpha of 0.3. You can specify `facecolor` and + `edgecolor` separately for greater control. + facecolor : matplotlib color specification + Color for the polygon fill. + edgecolor : matplotlib color specification + Color for the polygon boundary. + linewidth : float + The line width for the polygon boundary. + **kwargs + Additional keyword arguments passed to the matplotlib Patch. + + Returns + ------- + Matplotlib artist (PathPatch), if `add_points` is false. + A tuple of Matplotlib artists (PathPatch, Line2D), if `add_points` is true. + """ + from matplotlib import colors + + if ax is None: + ax = _default_ax() + + if color is None: + color = "C0" + color = colors.to_rgba(color) + + if facecolor is None: + facecolor = list(color) + facecolor[-1] = 0.3 + facecolor = tuple(facecolor) + + if edgecolor is None: + edgecolor = color + + patch = patch_from_polygon( + polygon, + facecolor=facecolor, + edgecolor=edgecolor, + linewidth=linewidth, + **kwargs, + ) + ax.add_patch(patch) + ax.autoscale_view() + + if add_points: + line = plot_points(polygon, ax=ax, color=color) + return patch, line + + return patch + + +def plot_points(geom, ax=None, color=None, marker="o", **kwargs): + """ + Plot a Point/MultiPoint or the vertices of any other geometry type. + + Parameters + ---------- + geom : shapely.Geometry + Any shapely Geometry object, from which all vertices are extracted + and plotted. + ax : matplotlib Axes, default None + The axes on which to draw the plot. If not specified, will get the + current active axes or create a new figure. + color : matplotlib color specification + Color for the filled points. You can use `markeredgecolor` and + `markerfacecolor` to have different edge and fill colors. + marker : str, default "o" + The matplotlib marker for the points. + **kwargs + Additional keyword arguments passed to matplotlib `plot` (Line2D). + + Returns + ------- + Matplotlib artist (Line2D) + """ + if ax is None: + ax = _default_ax() + + coords = shapely.get_coordinates(geom) + (line,) = ax.plot( + coords[:, 0], + coords[:, 1], + linestyle="", + marker=marker, + color=color, + **kwargs, + ) + return line + + +# patch shapely plotting capabilities if not found +if not hasattr(shapely, "plotting"): + shapely.plotting = type("plotting", (), {}) + shapely.plotting.plot_polygon = plot_polygon + shapely.plotting.plot_points = plot_points + shapely.plotting.patch_from_polygon = patch_from_polygon + shapely.plotting._path_from_polygon = _path_from_polygon +# --- END backport of shapely plotting capabilities --- + + +def _solar_vector(zenith, azimuth): + """ + Calculate the solar vector in 3D space. Origins from the Sun to the + observer on Earth characterized by the zenith and azimuth angles. + + Parameters + ---------- + zenith : numeric + Solar zenith angle. In degrees [°]. + azimuth : numeric + Solar azimuth angle. Positive is clockwise from the North in the + horizontal plane. North=0°, East=90°, South=180°, West=270°. + In degrees [°]. + + Returns + ------- + numpy.ndarray, shape (3,) or (N, 3) + Unitary solar vector. ``N`=len(zenith)=len(azimuth)``. + + References + ---------- + .. [1] E. Lorenzo, L. Narvarte, and J. Muñoz, 'Tracking and + back-tracking', Progress in Photovoltaics: Research and Applications, + vol. 19, no. 6, pp. 747-753, 2011, :doi:`10.1002/pip.1085`. + .. [2] Kevin S. Anderson, Adam R. Jensen; Shaded fraction and + backtracking in single-axis trackers on rolling terrain. J. Renewable + Sustainable Energy 1 March 2024; 16 (2): 023504. + :doi:`10.1063/5.0202220`. + """ + # Eq. (2), [1]; with zenith instead of elevation; coordinate system of [2] + return np.array( + [ + sind(zenith) * sind(azimuth), + sind(zenith) * cosd(azimuth), + cosd(zenith), + ] + ) + + +def _plane_normal_vector(tilt, azimuth): + """ + Calculate the normal vector of a plane defined by a tilt and azimuth. + + See Eq. (18) of [1]. It has been changed to match system coordinates in + Fig. 1 of [2]. + + Parameters + ---------- + azimuth : numeric + Azimuth angle of the plane. Positive is clockwise from the North in + the horizontal plane. North=0°, East=90°, South=180°, West=270°. + In degrees [°]. + tilt : numeric + Tilt angle of the plane. Positive is downwards from the horizontal in + the direction of ``azimuth``. In degrees [°]. + + Returns + ------- + numpy.ndarray, shape (3,) or (N, 3) + Unitary normal vector of the plane. ``N`=len(azimuth)=len(tilt)``. + + References + ---------- + .. [1] S. Zainali et al., 'Direct and diffuse shading factors modelling + for the most representative agrivoltaic system layouts', Applied + Energy, vol. 339, p. 120981, Jun. 2023, + :doi:`10.1016/j.apenergy.2023.120981`. + .. [2] Kevin S. Anderson, Adam R. Jensen; Shaded fraction and + backtracking in single-axis trackers on rolling terrain. J. Renewable + Sustainable Energy 1 March 2024; 16 (2): 023504. + :doi:`10.1063/5.0202220`. + """ + # Eq. (18) of [1], but coordinate system specified in Fig. 1, [2] + return np.array( + [sind(tilt) * sind(azimuth), sind(tilt) * cosd(azimuth), cosd(tilt)] + ) + + +class FlatSurface: + """ + Represents a flat surface in 3D space with a given azimuth and tilt and + boundaries defined by a shapely Polygon. + Allows to calculate the shading on this surface from other objects, + both in 2D and 3D. + In addition, it can + + .. warning:: + + This constructor does **not** check the ``azimuth`` and ``tilt`` match + the ``polygon`` orientation nor the ``polygon`` vertices are coplanar. + It is the user's responsibility to ensure the surface is correctly + defined. + + Parameters + ---------- + tilt : float + Surface tilt, angle it is inclined with respect to the horizontal + plane. Tilted downwards ``azimuth``. + 0°=Horizontal, 90°=Vertical. In degrees [°]. + azimuth : float + Surface azimuth, angle at which it points downwards. + 0°=North, 90°=East, 180°=South, 270°=West. In degrees [°]. + polygon : shapely.Polygon or array[N, 3] + Shapely Polygon or boundaries to build it. + Holes are ignored for now. + roll : float, default 0.0° + Right-handed rotation around the surface normal vector defined by + ``tilt`` and ``azimuth``. In degrees [°]. + reference_point : array-like, shape (3,), optional + Point to use as reference for 2D projections. If not provided, the + first vertex of the polygon is used. + + References + ---------- + .. [1] S. Zainali et al., 'Direct and diffuse shading factors modelling + for the most representative agrivoltaic system layouts', Applied + Energy, vol. 339, p. 120981, Jun. 2023, + :doi:`10.1016/j.apenergy.2023.120981`. + .. [2] Y. Cascone, V. Corrado, and V. Serra, 'Calculation procedure of + the shading factor under complex boundary conditions', Solar Energy, + vol. 85, no. 10, pp. 2524-2539, Oct. 2011, + :doi:`10.1016/j.solener.2011.07.011`. + .. [3] Kevin S. Anderson, Adam R. Jensen; Shaded fraction and + backtracking in single-axis trackers on rolling terrain. J. Renewable + Sustainable Energy 1 March 2024; 16 (2): 023504. + :doi:`10.1063/5.0202220`. + """ + + @property + def azimuth(self): + return self._azimuth + + @property + def tilt(self): + return self._tilt + + @property + def roll(self): + return self._roll + + @property + def polygon(self): + return self._polygon + + @property + def reference_point(self): + return self._reference_point + + def __init__( + self, + azimuth, + tilt, + polygon_boundaries, + roll=0.0, + reference_point=None, + ): + self._azimuth = np.array(azimuth) + self._tilt = np.array(tilt) + self._roll = np.array(roll) + # works for polygon_boundaries := array[N, 3] | shapely.Polygon + self._polygon = Polygon(polygon_boundaries) + # internal 2D coordinates-system to translate projections matrix + # only defined if needed later on + self._reference_point = ( + reference_point + if reference_point is not None + else self._polygon.exterior.coords[0] + ) + self._rotation_inverse = None + self._projected_polygon = None + + def __repr__(self): + return None + + def _transform_to_2D(self, points): + """ + Transform a 3D point that **belongs** to the surface, + to the 2D plane of this surface with respect to the local reference. + + Essentially, it undoes the surface rotations to make the third + coordinate zero. + + Parameters + ---------- + point : array-like, shape (N, 3) + Point to project. + + Returns + ------- + array-like, shape (N, 2) + Projected point. + + Raises + ------ + RuntimeError + If the 2D projection fails to transform to the surface plane. + Check the input points belong to the surface. + + References + ---------- + .. [1] Y. Cascone, V. Corrado, and V. Serra, 'Calculation procedure of + the shading factor under complex boundary conditions', Solar Energy, + vol. 85, no. 10, pp. 2524-2539, Oct. 2011, + :doi:`10.1016/j.solener.2011.07.011`. + """ + # undo surface rotations to make the third coordinate zero + _projection = self._get_projection_transform() + # Section 4.3 in [1], with a plane point + vertices_2d = _projection.apply(points - self.reference_point) + if not np.allclose(vertices_2d[:, 2], 0.0, atol=1e-10): + raise RuntimeError( + "Non-null third coordinate in 2D projection. " + + "This should not happen. Please report a bug." + ) + return np.delete(vertices_2d, 2, axis=1) + + def get_3D_shades_from(self, solar_zenith, solar_azimuth, *others): + """ + Calculate 3D shades on this surface from other objects. + + 3D shade points are guaranteed to be on the projected plane of this + surface. This method is useful to plot the resulting shade. + For the shades referred to the surface and a valid area property, + use the faster method :py:method:`get_2D_shades_from`. + + Shade is clipped to this object boundaries. + + Parameters + ---------- + solar_zenith : float + Solar zenith angle. In degrees [°]. + solar_azimuth : float + Solar azimuth angle. In degrees [°]. + others : FlatSurface or derived + Obstacles whose shadow will be projected onto this surface. + + Returns + ------- + shapely.MultiPolygon + Shapely Polygon objects representing the shades on this surface. + These shades may overlap. + + Notes + ----- + This method is based on the algorithm presented in [1]_ to calculate + the shades, but adds an extra step (a point that belongs to the plane) + to project the objects onto the surface plane instead of the plane that + goes through the origin. + + References + ---------- + .. [1] Y. Cascone, V. Corrado, and V. Serra, 'Calculation procedure of + the shading factor under complex boundary conditions', Solar Energy, + vol. 85, no. 10, pp. 2524-2539, Oct. 2011, + :doi:`10.1016/j.solener.2011.07.011`. + """ + solar_vec = solar_vec = _solar_vector( # Eq. (8) -> x,y,z + solar_zenith, solar_azimuth + ) + normal_vec = _plane_normal_vector( # Eq. (18) -> a,b,c + self._tilt, self._azimuth + ) + + def project_point_to_real_plane(vertex): # vertex -> Px, Py, Pz + # Similar to Eq. (20), but takes into account plane position so + # result belongs to the plane that contains the bounded surface + t = ((self._reference_point - vertex) @ normal_vec) / ( + solar_vec @ normal_vec + ) + p_prime = vertex + (t * solar_vec.T).T # Eq. (19) + return p_prime + + _polygon = self._polygon # intersect with 3D surface + + def get_3D_shade_from_flat_surface(other): + coords_to_project = np.array(other.polygon.exterior.coords[:-1]) + projected_vertices = np.fromiter( + map(project_point_to_real_plane, coords_to_project), + dtype=(float, 3), + count=len(coords_to_project), # speeds up allocation + ) + # create shapely shade object and bound it to the surface + try: # shapely > 2.0 + shade = Polygon(projected_vertices).intersection( + _polygon, grid_size=1e-12 + ) + except TypeError: # shapely < 2.0 + shade = Polygon(projected_vertices).intersection(_polygon) + return shade if isinstance(shade, Polygon) else Polygon() + + return MultiPolygon(map(get_3D_shade_from_flat_surface, others)) + + def get_2D_shades_from( + self, solar_zenith, solar_azimuth, *others, shades_3d=None + ): + """ + Calculate 2D shades on this surface from obstacles ``others*``. + + If ``shades`` is provided, the method will combine those shades with + the ones casted from the objects in ``others`` with the existing ones. + Use this parameter if you already have the 3D shades to avoid + recalculating them. + + Parameters + ---------- + solar_zenith : float + Solar zenith angle. In degrees [°]. + solar_azimuth : float + Solar azimuth angle. N=0°, E=90°, S=180°, W=270°. In degrees [°]. + others : FlatSurface or derived + Obstacles whose shadow will be projected onto this surface. + shades_3d : shapely.MultiPolygon, optional + Already calculated 3D shades to combine with the new ones, if any. + They can be obtained from :py:method:`get_3D_shades_from`. + + Returns + ------- + shapely.MultiPolygon + Collection of polygons representing the 2D shades clipped to + this surface. + + Raises + ------ + RuntimeError + If the 2D projection of the shades fails to transform to the + surface plane. In this case, you must ensure ``shades_3d`` is + correctly calculated and those shades belong to the object you + are projecting onto. + + Notes + ----- + This method is based on the algorithm presented in [1]_ to calculate + the shades, but adds an extra step to translate the projection objects + to the plane that goes through the origin, to remove the Z coordinate + and make the projection in 2D space. + + References + ---------- + .. [1] Y. Cascone, V. Corrado, and V. Serra, 'Calculation procedure of + the shading factor under complex boundary conditions', Solar Energy, + vol. 85, no. 10, pp. 2524-2539, Oct. 2011, + :doi:`10.1016/j.solener.2011.07.011`. + """ + shades = self.get_3D_shades_from(solar_zenith, solar_azimuth, *others) + if shades_3d is not None: + shades = MultiPolygon((*shades.geoms, *shades_3d.geoms)) + + # and clip the 2D shades to the surface in 2D + _self_projected_polygon = self.representation_in_2D_space() + + def get_2D_shade_from_flat_surface(other): + coords_to_project = np.array(other.exterior.coords[:-1]) + + # create shapely shade object and bound it to the surface + projected_vertices_2d = self._transform_to_2D(coords_to_project) + shade = Polygon(projected_vertices_2d).intersection( + _self_projected_polygon + ) + # Return Polygons only, geometries with empty area are omitted + return shade if isinstance(shade, Polygon) else Polygon() + + return MultiPolygon(map(get_2D_shade_from_flat_surface, shades.geoms)) + + def _get_projection_transform(self): + """ + Get the transformation matrix to project points to the 2D plane of + this surface. + + Returns + ------- + scipy.spatial.transform.Rotation + Transformation matrix to project points to the 2D plane of this + surface. + + References + ---------- + .. [1] Y. Cascone, V. Corrado, and V. Serra, 'Calculation procedure of + the shading factor under complex boundary conditions', Solar Energy, + vol. 85, no. 10, pp. 2524-2539, Oct. 2011, + :doi:`10.1016/j.solener.2011.07.011`. + """ + # Section 4.3 in [1], plus another rotation to align the surface: + # the rotation around the normal vector + if self._rotation_inverse is None: + self._rotation_inverse = Rotation.from_euler( + "zxz", + [self._azimuth - 180, -self._tilt, -self._roll], + degrees=True, + ) + return self._rotation_inverse + + def representation_in_2D_space(self): + """ + Get the 2D representation of the surface in its local plane, with + respect to the reference point. + + Returns + ------- + shapely.Polygon + 2D representation of the surface. + + Notes + ----- + This method is useful to plot the surface in 2D space and see + where the shadows are casted. + """ + if self._projected_polygon is None: + vertices = np.array(self._polygon.exterior.coords[:-1]) + self._projected_polygon = Polygon(self._transform_to_2D(vertices)) + return self._projected_polygon + + def combine_2D_shades(self, *shades): + """ + Combine overlapping shades into a single one, but keep non-overlapping + shades separated. + + Parameters + ---------- + shades : shapely.MultiPolygon + Shapely MultiPolygon object representing the shades. + + Returns + ------- + shapely.MultiPolygon + Combined shades. + """ + # TODO: implement this method + return None + + def plot(self, ax=None, **kwargs): + pass # TODO: implement this method + + +class RectangularSurface(FlatSurface): + """ + Represents a rectangular surface in 3D space with a given ``azimuth``, + ``tilt``, ``axis_tilt`` and a center point. This is a subclass of + :py:class:`~pvlib.spatial.FlatSurface` handy for rectangular surfaces + like PV arrays. + + See :py:class:`pvlib.spatial.FlatSurface` for information on methods and + properties. + + Parameters + ---------- + center : array-like, shape (3,) + Center of the surface + azimuth : float + Azimuth of the surface. Positive is clockwise from the North in the + horizontal plane. North=0°, East=90°, South=180°, West=270°. + In degrees [°]. + tilt : float + Tilt of the surface, angle it is inclined with respect to the + horizontal plane. Positive is downwards ``azimuth``. + In degrees [°]. + axis_tilt : float + Rotation around the normal vector of the surface. Right-handed. + In degrees [°]. + width: width of the surface + For a horizontal surface, the width is parallel to the azimuth + length: length of the surface + Perpendicular to the surface azimuth + reference_point : array-like, shape (3,), optional + Point to use as reference for 2D projections. If not provided, the + central point is used. + """ + + def __init__( + self, + center, + azimuth, + tilt, + axis_tilt, + width, + length, + reference_point=None, + ): + corners = np.array( + [ + [-length / 2, -width / 2, 0], + [-length / 2, +width / 2, 0], + [+length / 2, +width / 2, 0], + [+length / 2, -width / 2, 0], + ] + ) + # rotate corners to match the surface orientation + # note pvlib convention uses a left-handed azimuth rotation + # and tilt is defined as the angle from the horizontal plane + # downwards the azimuth direction (also left-handed) + # axis_tilt is the rotation around the normal vector, right-handed + _rotation = Rotation.from_euler( + "ZXZ", [-azimuth, -tilt, axis_tilt], degrees=True + ) + polygon = Polygon(_rotation.apply(corners) + center) + super().__init__( + azimuth=azimuth, + tilt=tilt, + polygon_boundaries=polygon, + roll=axis_tilt, + reference_point=reference_point if reference_point else center, + ) + + def plot(self, ax=None, **kwargs): + """ + Plot the rectangular surface. + + Parameters + ---------- + ax : mpl_toolkits.mplot3d.axes3d.Axes3D, optional + Axes where to plot the surface. If None, a new figure is created. + **kwargs : dict + Additional arguments passed to + :external:ref:`mpl_toolkits.mplot3d.axes3d.Axes3D.plot_trisurf`. + """ + if ax is None: + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + x, y, z = np.hsplit( + np.array(self._polygon.exterior.coords[:-1]).flatten(order="F"), 3 + ) + return ax.plot_trisurf( + x, y, z, triangles=((0, 1, 2), (0, 2, 3)), **kwargs + ) diff --git a/pvlib/tests/test_spatial.py b/pvlib/tests/test_spatial.py new file mode 100644 index 0000000000..7ddcd5eed0 --- /dev/null +++ b/pvlib/tests/test_spatial.py @@ -0,0 +1,174 @@ +"""Tests spatial submodule.""" + +from pvlib import spatial +from numpy.testing import assert_allclose +from shapely.geometry import Polygon + +import pytest + + +@pytest.mark.parametrize( + "azimuth, zenith, expected", + [ + # vertical vector + (0, 0, [0, 0, 1]), + (90, 0, [0, 0, 1]), + (180, 0, [0, 0, 1]), + (270, 0, [0, 0, 1]), + (0, 90, [0, 1, 0]), + (90, 90, [1, 0, 0]), + (180, 90, [0, -1, 0]), + (270, 90, [-1, 0, 0]), + (0, 45, [0, 0.70710678, 0.70710678]), + (90, 45, [0.70710678, 0, 0.70710678]), + (180, 45, [0, -0.70710678, 0.70710678]), + (270, 45, [-0.70710678, 0, 0.70710678]), + ], +) +def test__solar_vector(azimuth, zenith, expected): + vec = spatial._solar_vector(zenith=zenith, azimuth=azimuth) + assert_allclose(vec, expected, atol=1e-15) + + +@pytest.mark.parametrize( + "azimuth, tilt, expected", + [ + # horizontal surfaces + (0, 0, [0, 0, 1]), + (90, 0, [0, 0, 1]), + (180, 0, [0, 0, 1]), + (270, 0, [0, 0, 1]), + # vertical surfaces + (0, 90, [0, 1, 0]), + (90, 90, [1, 0, 0]), + (180, 90, [0, -1, 0]), + (270, 90, [-1, 0, 0]), + # tilted surfaces + (0, 45, [0, 0.70710678, 0.70710678]), + (90, 45, [0.70710678, 0, 0.70710678]), + (180, 45, [0, -0.70710678, 0.70710678]), + (270, 45, [-0.70710678, 0, 0.70710678]), + ], +) +def test__plane_normal_vector(azimuth, tilt, expected): + vec = spatial._plane_normal_vector(tilt=tilt, azimuth=azimuth) + assert_allclose(vec, expected, atol=1e-8) + + +def test_FlatSurface__init__(): + # construct with native shapely polygon + surface_azimuth = 180 + surface_tilt = 0 + polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) + surface = spatial.FlatSurface( + azimuth=surface_azimuth, + tilt=surface_tilt, + polygon_boundaries=polygon, + ) + assert surface.azimuth == surface_azimuth + assert surface.tilt == surface_tilt + assert surface.polygon == polygon + assert isinstance(surface.polygon, Polygon) + # construct from coordinates + surface_azimuth = 180 + surface_tilt = 0 + polygon = [(0, 0), (1, 0), (1, 1), (0, 1)] + surface = spatial.FlatSurface( + azimuth=surface_azimuth, tilt=surface_tilt, polygon_boundaries=polygon + ) + assert surface.azimuth == surface_azimuth + assert surface.tilt == surface_tilt + assert surface.polygon == Polygon(polygon) + assert isinstance(surface.polygon, Polygon) + + +def test_FlatSurface_readonly_properties(): + surface_azimuth = 180 + surface_tilt = 0 + polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) + surface = spatial.FlatSurface( + azimuth=surface_azimuth, + tilt=surface_tilt, + polygon_boundaries=polygon, + ) + with pytest.raises(AttributeError): + surface.azimuth = 0 + with pytest.raises(AttributeError): + surface.tilt = 0 + with pytest.raises(AttributeError): + surface.polygon = polygon + + +def test_FlatSurface_get_3D_shades_from(): + pass + + +def test_FlatSurface_get_2D_shades_from(): + pass + + +def test_FlatSurface_combine_2D_shades(): + pass + + +def test_FlatSurface_plot(): + pass + + +def test_RectangularSurface__init__(): + # construct with native shapely polygon + center = [0, 0, 0] + surface_azimuth = 180 + surface_tilt = 0 + axis_tilt = 0 + width = 1 + length = 1 + surface = spatial.RectangularSurface( + center=center, + azimuth=surface_azimuth, + tilt=surface_tilt, + axis_tilt=axis_tilt, + width=width, + length=length, + ) + assert surface.reference_point == center + assert surface.azimuth == surface_azimuth + assert surface.tilt == surface_tilt + assert surface.roll == axis_tilt + # construct from coordinates + center = [0, 0, 0] + surface_azimuth = 180 + surface_tilt = 0 + axis_tilt = 0 + width = 1 + length = 1 + surface = spatial.RectangularSurface( + center=center, + azimuth=surface_azimuth, + tilt=surface_tilt, + axis_tilt=axis_tilt, + width=width, + length=length, + ) + assert surface.reference_point == center + assert surface.azimuth == surface_azimuth + assert surface.tilt == surface_tilt + assert surface.roll == axis_tilt + + +def test_RectangularSurface__calc_surface_tilt_and_azimuth(): + center = [0, 0, 0] + width = length = 2 + azimuth = 180 + tilt = 45 + axis_tilt = 45 + surface = spatial.RectangularSurface( + center=center, + azimuth=azimuth, + tilt=tilt, + axis_tilt=axis_tilt, + width=width, + length=length, + ) + assert surface.tilt == tilt + assert surface.azimuth == azimuth diff --git a/pvlib/tools.py b/pvlib/tools.py index 3d766b6f72..e26d158ff2 100644 --- a/pvlib/tools.py +++ b/pvlib/tools.py @@ -117,6 +117,23 @@ def atand(number): return res +def atan2d(y, x): + """ + Trigonometric quadrant-wise inverse tangent returning an angle in degrees. + + Parameters + ---------- + y, x : numeric + Input numbers + + Returns + ------- + result : numeric + arctan2 result in degrees + """ + return np.degrees(np.arctan2(y, x)) + + def localize_to_utc(time, location): """ Converts or localizes a time series to UTC. diff --git a/pyproject.toml b/pyproject.toml index de51b7c22c..a96d2883ad 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,6 +53,7 @@ optional = [ 'numba >= 0.17.0', 'solarfactors', 'statsmodels', + 'shapely >= 2.0.2', ] doc = [ 'ipython', @@ -64,6 +65,7 @@ doc = [ 'pillow', 'sphinx-toggleprompt >= 0.0.5', 'solarfactors', + 'shapely', ] test = [ 'pytest',