From 46a6de95812ac2d578cf2d3433d45bfb27d16d7e Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Mon, 6 May 2024 00:18:54 +0200 Subject: [PATCH 01/25] par shading factor dev commit - to be reverted --- .gitignore | 1 + PAR_direct_and_diffuse_sfs.py | 148 ++++++++++++++++++++++++++++++++++ 2 files changed, 149 insertions(+) create mode 100644 PAR_direct_and_diffuse_sfs.py diff --git a/.gitignore b/.gitignore index e124a08fdb..6465a1d3bb 100644 --- a/.gitignore +++ b/.gitignore @@ -98,3 +98,4 @@ coverage.xml env results +/.venv diff --git a/PAR_direct_and_diffuse_sfs.py b/PAR_direct_and_diffuse_sfs.py new file mode 100644 index 0000000000..ae1f3d020b --- /dev/null +++ b/PAR_direct_and_diffuse_sfs.py @@ -0,0 +1,148 @@ +# %% +import matplotlib.pyplot as plt +import numpy as np + +# %% +x0, y0, z0 = 0, 0, 1 / 2 # m +W, L = 1, 2 # m +corners0 = np.array( + [ + [x0, y0, z0], + [x0, y0 + W, z0], + [x0 + L, y0, z0], + [x0 + L, y0 + W, z0], + ] +) + + +# %% +# def rotation +def Ry(theta): + return np.array( + [ + [np.cos(theta), 0, np.sin(theta)], + [0, 1, 0], + [-np.sin(theta), 0, np.cos(theta)], + ] + ) + + +def Rx(theta): + return np.array( + [ + [1, 0, 0], + [0, np.cos(theta), -np.sin(theta)], + [0, np.sin(theta), np.cos(theta)], + ] + ) + + +# %% +# rotate corners, general abstraction of row +def rotate_row_with_center_pivot(corners, theta_x, theta_y): + Rot_x, Rot_y = Rx(theta_x), Ry(theta_y) + midpoint = np.mean(corners, axis=0) + Rot = Rot_x @ Rot_y + return np.array([Rot @ (corner - midpoint) + midpoint for corner in corners]) + + +# %% +# unittest previous function +def test_rotate_row_with_center_pivot(): + corners = np.array( + [ + [0, 0, 0], + [0, 1, 0], + [1, 0, 0], + [1, 1, 0], + ] + ) + corners_rotated = rotate_row_with_center_pivot(corners, -np.pi / 2, 0) + print(corners_rotated) + + +test_rotate_row_with_center_pivot() + + +# %% +def plot_corners(corners, fig=None, ax=None, **kwargs) -> plt.Figure: + if fig is None: + fig = plt.figure() + if ax is None: + ax = fig.add_subplot(111, projection="3d") + x_ = corners[:, 0].reshape(2, 2) + y_ = corners[:, 1].reshape(2, 2) + z_ = corners[:, 2].reshape(2, 2) + ax.plot_surface(x_, y_, z_, **kwargs) + return fig + + +# %% +fig = plot_corners(corners0, alpha=0.5) +fig.show() + +# %% +rotated_fixed_tilt = rotate_row_with_center_pivot(corners0, -np.pi / 2, 0) +fig = plot_corners(rotated_fixed_tilt, alpha=0.5) +fig.show() + + +# %% +def solar_vector(zenith, azimuth): + # Eq. (8), but in terms of zenith instead of elevation + return np.array( + [ + np.sin(zenith) * np.sin(azimuth), + np.sin(zenith) * np.cos(azimuth), + np.cos(zenith), + ] + ) + + +def normal_vector_of_row(tilt, azimuth): + # Eq. (18) + return np.array( + [ + np.sin(tilt) * np.cos(azimuth), + np.sin(tilt) * np.sin(azimuth), + np.cos(tilt), + ] + ) + + +def corners_projection_onto_plane( + corners, solar_zenith, solar_azimuth, plane_tilt, plane_azimuth +): + x, y, z = solar_vector(solar_zenith, solar_azimuth) + a, b, c = normal_vector_of_row(plane_tilt, plane_azimuth) + def projection(corner): + Px, Py, Pz = corner + # Eq. (20) + t = - (a*Px+b*Py+c*Pz) / (a*x+b*y+c*z) + # Eq. (19) + p_prime = np.array([Px + x*t, Py + y*t, Pz + z*t]) + return p_prime + return np.array([projection(corner) for corner in corners]) + + +# %% +# set hypothetical values +solar_zenith = np.pi / 2 - 1e-1 +solar_azimuth = np.pi +row_tilt = np.pi / 6 +row_azimuth = np.pi + +# create the row +corners1 = rotate_row_with_center_pivot(corners0, row_tilt, 0) +corners_projected = corners_projection_onto_plane(corners1, solar_zenith, solar_azimuth, 0, 0) + +fig = plt.figure() +ax = fig.add_subplot(111, projection="3d") +ax.set_xlim(0, 2) +ax.set_ylim(0, 3) +ax.set_zlim(0, 1) +_ = plot_corners(corners1, ax=ax, alpha=0.5, color="grey") +_ = plot_corners(corners_projected, ax=ax, alpha=0.5, color="brown") +fig.show() + +# %% From e483e92df952fc922b387204d17f0dd05f21bdd3 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sun, 26 May 2024 13:19:39 +0200 Subject: [PATCH 02/25] par shading factor dev commit - to revert --- PAR_direct_and_diffuse_sfs.py | 75 +++++++++++++++++++++++++---------- 1 file changed, 55 insertions(+), 20 deletions(-) diff --git a/PAR_direct_and_diffuse_sfs.py b/PAR_direct_and_diffuse_sfs.py index ae1f3d020b..293535dde6 100644 --- a/PAR_direct_and_diffuse_sfs.py +++ b/PAR_direct_and_diffuse_sfs.py @@ -1,6 +1,11 @@ # %% import matplotlib.pyplot as plt import numpy as np +from scipy.spatial import ConvexHull +from scipy.spatial.transform import Rotation as R +from pvlib.tools import sind, cosd + +from itertools import repeat # %% x0, y0, z0 = 0, 0, 1 / 2 # m @@ -17,12 +22,22 @@ # %% # def rotation +def Rz(theta): + return np.array( + [ + [cosd(theta), -sind(theta), 0], + [sind(theta), cosd(theta), 0], + [0, 0, 1], + ] + ) + + def Ry(theta): return np.array( [ - [np.cos(theta), 0, np.sin(theta)], + [cosd(theta), 0, sind(theta)], [0, 1, 0], - [-np.sin(theta), 0, np.cos(theta)], + [-sind(theta), 0, cosd(theta)], ] ) @@ -31,8 +46,8 @@ def Rx(theta): return np.array( [ [1, 0, 0], - [0, np.cos(theta), -np.sin(theta)], - [0, np.sin(theta), np.cos(theta)], + [0, cosd(theta), -sind(theta)], + [0, sind(theta), cosd(theta)], ] ) @@ -57,7 +72,7 @@ def test_rotate_row_with_center_pivot(): [1, 1, 0], ] ) - corners_rotated = rotate_row_with_center_pivot(corners, -np.pi / 2, 0) + corners_rotated = rotate_row_with_center_pivot(corners, -90, 0) print(corners_rotated) @@ -70,9 +85,9 @@ def plot_corners(corners, fig=None, ax=None, **kwargs) -> plt.Figure: fig = plt.figure() if ax is None: ax = fig.add_subplot(111, projection="3d") - x_ = corners[:, 0].reshape(2, 2) - y_ = corners[:, 1].reshape(2, 2) - z_ = corners[:, 2].reshape(2, 2) + x_ = corners[:, 0].reshape(2, -1) + y_ = corners[:, 1].reshape(2, -1) + z_ = corners[:, 2].reshape(2, -1) ax.plot_surface(x_, y_, z_, **kwargs) return fig @@ -82,7 +97,7 @@ def plot_corners(corners, fig=None, ax=None, **kwargs) -> plt.Figure: fig.show() # %% -rotated_fixed_tilt = rotate_row_with_center_pivot(corners0, -np.pi / 2, 0) +rotated_fixed_tilt = rotate_row_with_center_pivot(corners0, -90, 0) fig = plot_corners(rotated_fixed_tilt, alpha=0.5) fig.show() @@ -92,9 +107,9 @@ def solar_vector(zenith, azimuth): # Eq. (8), but in terms of zenith instead of elevation return np.array( [ - np.sin(zenith) * np.sin(azimuth), - np.sin(zenith) * np.cos(azimuth), - np.cos(zenith), + sind(zenith) * sind(azimuth), + sind(zenith) * cosd(azimuth), + cosd(zenith), ] ) @@ -103,9 +118,9 @@ def normal_vector_of_row(tilt, azimuth): # Eq. (18) return np.array( [ - np.sin(tilt) * np.cos(azimuth), - np.sin(tilt) * np.sin(azimuth), - np.cos(tilt), + sind(tilt) * cosd(azimuth), + sind(tilt) * sind(azimuth), + cosd(tilt), ] ) @@ -127,14 +142,16 @@ def projection(corner): # %% # set hypothetical values -solar_zenith = np.pi / 2 - 1e-1 -solar_azimuth = np.pi -row_tilt = np.pi / 6 -row_azimuth = np.pi +solar_zenith = 30 +solar_azimuth = 180 +row_tilt = 30 +row_azimuth = 180 +plane_tilt = -10 +plane_azimuth = 90 # create the row corners1 = rotate_row_with_center_pivot(corners0, row_tilt, 0) -corners_projected = corners_projection_onto_plane(corners1, solar_zenith, solar_azimuth, 0, 0) +corners_projected = corners_projection_onto_plane(corners1, solar_zenith, solar_azimuth, plane_tilt, plane_azimuth) fig = plt.figure() ax = fig.add_subplot(111, projection="3d") @@ -146,3 +163,21 @@ def projection(corner): fig.show() # %% +# Compare corners_projected with corners after rotation to new coordinate systems of the projection plane, to remove the third coordinate + +def from_3d_plane_to_2d(points, tilt, azimuth): + # Section 4.3 in [2] + R_z = Rz(90 + azimuth) + R_x = Rx(tilt) + rot = (R_z @ R_x).T + if points.ndim == 1: + points = points.reshape(1, -1) + new_points = np.fromiter(map(rot.__matmul__, points), dtype=(float, 3)) + return np.delete(new_points, 2, axis=1) + + +new_points = from_3d_plane_to_2d(corners_projected, plane_tilt, plane_azimuth) + +print(new_points) + +# %% From 1dc536622c7e763650349baa08b3bc46cafb754c Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sun, 26 May 2024 23:23:22 +0200 Subject: [PATCH 03/25] Update PAR_direct_and_diffuse_sfs.py --- PAR_direct_and_diffuse_sfs.py | 1 + 1 file changed, 1 insertion(+) diff --git a/PAR_direct_and_diffuse_sfs.py b/PAR_direct_and_diffuse_sfs.py index 293535dde6..5a76fa3872 100644 --- a/PAR_direct_and_diffuse_sfs.py +++ b/PAR_direct_and_diffuse_sfs.py @@ -173,6 +173,7 @@ def from_3d_plane_to_2d(points, tilt, azimuth): if points.ndim == 1: points = points.reshape(1, -1) new_points = np.fromiter(map(rot.__matmul__, points), dtype=(float, 3)) + assert np.allclose(new_points[:, 2], 0), "The third coordinate should be zero, check input parameters are consistent with the projection plane." return np.delete(new_points, 2, axis=1) From 44fb93975ffa96e7ccc576203b30efdc488e1b39 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Thu, 30 May 2024 23:26:29 +0200 Subject: [PATCH 04/25] Make shapely an optional dependency --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index de51b7c22c..273cf2a2a7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,6 +53,7 @@ optional = [ 'numba >= 0.17.0', 'solarfactors', 'statsmodels', + 'shapely', ] doc = [ 'ipython', From a9ee106d11add183a68259d34b32cd1e92ec2415 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Mon, 3 Jun 2024 01:49:13 +0200 Subject: [PATCH 05/25] Update PAR_direct_and_diffuse_sfs.py --- PAR_direct_and_diffuse_sfs.py | 195 +++++++++++++++++++++++++++++++++- 1 file changed, 190 insertions(+), 5 deletions(-) diff --git a/PAR_direct_and_diffuse_sfs.py b/PAR_direct_and_diffuse_sfs.py index 5a76fa3872..c06fef1872 100644 --- a/PAR_direct_and_diffuse_sfs.py +++ b/PAR_direct_and_diffuse_sfs.py @@ -58,7 +58,9 @@ def rotate_row_with_center_pivot(corners, theta_x, theta_y): Rot_x, Rot_y = Rx(theta_x), Ry(theta_y) midpoint = np.mean(corners, axis=0) Rot = Rot_x @ Rot_y - return np.array([Rot @ (corner - midpoint) + midpoint for corner in corners]) + return np.array( + [Rot @ (corner - midpoint) + midpoint for corner in corners] + ) # %% @@ -130,13 +132,15 @@ def corners_projection_onto_plane( ): x, y, z = solar_vector(solar_zenith, solar_azimuth) a, b, c = normal_vector_of_row(plane_tilt, plane_azimuth) + def projection(corner): Px, Py, Pz = corner # Eq. (20) - t = - (a*Px+b*Py+c*Pz) / (a*x+b*y+c*z) + t = -(a * Px + b * Py + c * Pz) / (a * x + b * y + c * z) # Eq. (19) - p_prime = np.array([Px + x*t, Py + y*t, Pz + z*t]) + p_prime = np.array([Px + x * t, Py + y * t, Pz + z * t]) return p_prime + return np.array([projection(corner) for corner in corners]) @@ -151,7 +155,9 @@ def projection(corner): # create the row corners1 = rotate_row_with_center_pivot(corners0, row_tilt, 0) -corners_projected = corners_projection_onto_plane(corners1, solar_zenith, solar_azimuth, plane_tilt, plane_azimuth) +corners_projected = corners_projection_onto_plane( + corners1, solar_zenith, solar_azimuth, plane_tilt, plane_azimuth +) fig = plt.figure() ax = fig.add_subplot(111, projection="3d") @@ -165,6 +171,7 @@ def projection(corner): # %% # Compare corners_projected with corners after rotation to new coordinate systems of the projection plane, to remove the third coordinate + def from_3d_plane_to_2d(points, tilt, azimuth): # Section 4.3 in [2] R_z = Rz(90 + azimuth) @@ -173,7 +180,9 @@ def from_3d_plane_to_2d(points, tilt, azimuth): if points.ndim == 1: points = points.reshape(1, -1) new_points = np.fromiter(map(rot.__matmul__, points), dtype=(float, 3)) - assert np.allclose(new_points[:, 2], 0), "The third coordinate should be zero, check input parameters are consistent with the projection plane." + assert np.allclose( + new_points[:, 2], 0 + ), "The third coordinate should be zero, check input parameters are consistent with the projection plane." return np.delete(new_points, 2, axis=1) @@ -181,4 +190,180 @@ def from_3d_plane_to_2d(points, tilt, azimuth): print(new_points) +# %% +########################## +# OOP Implementation +########################## +import numpy as np +import shapely as sp +from pvlib.tools import sind, cosd + +from abc import ABC, abstractmethod + + +def Rz(theta): + return np.array( + [ + [cosd(theta), -sind(theta), 0], + [sind(theta), cosd(theta), 0], + [0, 0, 1], + ] + ) + + +def Ry(theta): + return np.array( + [ + [cosd(theta), 0, sind(theta)], + [0, 1, 0], + [-sind(theta), 0, cosd(theta)], + ] + ) + + +def Rx(theta): + return np.array( + [ + [1, 0, 0], + [0, cosd(theta), -sind(theta)], + [0, sind(theta), cosd(theta)], + ] + ) + + +def solar_vector(zenith, azimuth): + # Eq. (8), but in terms of zenith instead of elevation + return np.array( + [ + sind(zenith) * sind(azimuth), + sind(zenith) * cosd(azimuth), + cosd(zenith), + ] + ) + + +def normal_vector_of_surface(tilt, azimuth): + # Eq. (18) + return np.array( + [ + sind(tilt) * cosd(azimuth), + sind(tilt) * sind(azimuth), + cosd(tilt), + ] + ) + + +class CoplanarSurface: + def __init__(self, azimuth, tilt, polygon_boundaries, point_on_surface=None): + """ + + .. warning:: + + This constructor does not make any check on the input parameters. + Use any of the class methods ``from_*`` to create an object. + + Parameters + ---------- + azimuth : float + Surface azimuth, angle at which it points downwards. 0° is North, 90° is East, 180° is South, 270° is West. In degrees [°]. + tilt : float + Surface tilt, angle at which it is inclined with respect to the horizontal plane. Tilted downwards ``azimuth``. 0° is horizontal, 90° is vertical. In degrees [°]. + polygon : shapely.LinearRing or array[N, 3] + Shapely Polygon or boundaries to build it. If you need to specify holes, provide your custom LinearRing. Note that not all shapely objects calculate a ``point_on_surface`` correctly in 3D space. LinearRing is guaranteed to work. + """ + self.azimuth = azimuth + self.tilt = tilt + # works for polygon_boundaries := array[N, 3] | shapely.Polygon + self.polygon = sp.LinearRing(polygon_boundaries) + # representative point to calculate if obstacles are in front or behind + if point_on_surface is None: + self.belonging_point = self.polygon.point_on_surface() + else: + self.belonging_point = point_on_surface + # internal 2D coordinates-system to translate projections + # TODO + + + def get_shade_by(self, solar_zenith, solar_azimuth, *others): + # project ``others`` vertices onto this shaded surface + solar_vec = solar_vector(solar_zenith, solar_azimuth) # := x, y, z + normal_vec = normal_vector_of_surface(self.tilt, self.azimuth) # := a, b, c + def projection(corner): # corner := Px, Py, Pz + # Eq. (20) + t = -(normal_vec * corner) / (solar_vec * normal_vec) + # Eq. (19) + p_prime = corner + t * solar_vec # Px + x * t, Py + y * t, Pz + z * t + return p_prime + corns = np.array([projection(corner) for corner in self.polygon.coords[:-1]]) + # undo surface rotations to make the third coordinate zero + # TODO + new_points = from_3d_plane_to_2d( + corners_projected, + self.base.surface_tilt, + self.base.surface_azimuth, + ) + # create a polygon + polygon = sp.geometry.Polygon(new_points) + # create a convex hull + hull = ConvexHull(new_points) + # create a shapely polygon + shapely_hull = sp.geometry.Polygon(hull.points[hull.vertices]) + return polygon, shapely_hull + + +class RectangularSurface(CoplanarSurface): + def __init__( + self, center, surface_azimuth, surface_tilt, axis_tilt, width, length + ): + """ + + .. warning:: + + This constructor does not make any check on the input parameters. + Use any of the class methods ``from_*`` to create an object. + + center: center of the surface + surface_azimuth: azimuth of the surface + surface_tilt: tilt of the surface + width: width of the surface + length: length of the surface + TODO: which one is perpendicular to azimuth? + """ + self.center = np.array(center) + x_c, y_c, z_c = center + self.surface_azimuth = surface_azimuth + self.surface_tilt = surface_tilt + self.axis_tilt = axis_tilt + self.width = width + self.length = length + corners = np.array( + [ + [x_c - length / 2, y_c - width / 2, z_c], + [x_c - length / 2, y_c + width / 2, z_c], + [x_c + length / 2, y_c + width / 2, z_c], + [x_c + length / 2, y_c - width / 2, z_c], + ] + ) + # rotate corners to match the surface orientation + rX, rY, rZ = Rx(row_tilt), Ry(row_azimuth), Rz(axis_tilt) + rot = rX @ rY @ rZ + self.corners = np.array( + [rot @ (corner - self.center) + self.center for corner in corners] + ) + self.shapely_obj = sp.LinearRing( + [rot @ (corner - self.center) + self.center for corner in corners] + ) + + # @classmethod + # def from_slope_values(cls, center, slope_azimuth, slope_tilt, width, length): + # return cls(center, slope_azimuth, slope_tilt, width, length) + + # @classmethod + # def from_Array_values(cls, center, surface_azimuth, surface_tilt, axis_tilt, width, length): + # return cls(center, surface_azimuth, surface_tilt, axis_tilt, width, length) + + + + + # %% From 4163bb3aa14d0b299b4523311c23c1b35f104eb4 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Fri, 7 Jun 2024 16:56:16 +0200 Subject: [PATCH 06/25] Add atan2d --- pvlib/tools.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/pvlib/tools.py b/pvlib/tools.py index adf502a79d..e717e407e4 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. From 296153d5e3c55454b650a3179873698d278b6722 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Fri, 7 Jun 2024 16:56:19 +0200 Subject: [PATCH 07/25] Update PAR_direct_and_diffuse_sfs.py --- PAR_direct_and_diffuse_sfs.py | 226 ++++++++++++++++++++++++++-------- 1 file changed, 176 insertions(+), 50 deletions(-) diff --git a/PAR_direct_and_diffuse_sfs.py b/PAR_direct_and_diffuse_sfs.py index c06fef1872..e55b6d1b21 100644 --- a/PAR_direct_and_diffuse_sfs.py +++ b/PAR_direct_and_diffuse_sfs.py @@ -87,6 +87,8 @@ def plot_corners(corners, fig=None, ax=None, **kwargs) -> plt.Figure: fig = plt.figure() if ax is None: ax = fig.add_subplot(111, projection="3d") + if not isinstance(corners, np.ndarray): + corners = np.array(corners) x_ = corners[:, 0].reshape(2, -1) y_ = corners[:, 1].reshape(2, -1) z_ = corners[:, 2].reshape(2, -1) @@ -191,16 +193,23 @@ def from_3d_plane_to_2d(points, tilt, azimuth): print(new_points) # %% +# ########################## # OOP Implementation ########################## import numpy as np +import pandas as pd import shapely as sp -from pvlib.tools import sind, cosd +from pvlib.tools import sind, cosd, acosd +import matplotlib.pyplot as plt from abc import ABC, abstractmethod +def atan2d(y, x): + return np.degrees(np.arctan2(y, x)) + + def Rz(theta): return np.array( [ @@ -239,7 +248,7 @@ def solar_vector(zenith, azimuth): sind(zenith) * cosd(azimuth), cosd(zenith), ] - ) + ).T def normal_vector_of_surface(tilt, azimuth): @@ -250,11 +259,14 @@ def normal_vector_of_surface(tilt, azimuth): sind(tilt) * sind(azimuth), cosd(tilt), ] - ) + ).T -class CoplanarSurface: - def __init__(self, azimuth, tilt, polygon_boundaries, point_on_surface=None): +# %% +class FlatSurface: + def __init__( + self, azimuth, tilt, polygon_boundaries + ): """ .. warning:: @@ -269,49 +281,83 @@ def __init__(self, azimuth, tilt, polygon_boundaries, point_on_surface=None): tilt : float Surface tilt, angle at which it is inclined with respect to the horizontal plane. Tilted downwards ``azimuth``. 0° is horizontal, 90° is vertical. In degrees [°]. polygon : shapely.LinearRing or array[N, 3] - Shapely Polygon or boundaries to build it. If you need to specify holes, provide your custom LinearRing. Note that not all shapely objects calculate a ``point_on_surface`` correctly in 3D space. LinearRing is guaranteed to work. + Shapely Polygon or boundaries to build it. If you need to specify holes, provide your custom LinearRing. """ self.azimuth = azimuth self.tilt = tilt # works for polygon_boundaries := array[N, 3] | shapely.Polygon self.polygon = sp.LinearRing(polygon_boundaries) - # representative point to calculate if obstacles are in front or behind - if point_on_surface is None: - self.belonging_point = self.polygon.point_on_surface() - else: - self.belonging_point = point_on_surface - # internal 2D coordinates-system to translate projections - # TODO - - - def get_shade_by(self, solar_zenith, solar_azimuth, *others): - # project ``others`` vertices onto this shaded surface - solar_vec = solar_vector(solar_zenith, solar_azimuth) # := x, y, z - normal_vec = normal_vector_of_surface(self.tilt, self.azimuth) # := a, b, c - def projection(corner): # corner := Px, Py, Pz - # Eq. (20) - t = -(normal_vec * corner) / (solar_vec * normal_vec) - # Eq. (19) - p_prime = corner + t * solar_vec # Px + x * t, Py + y * t, Pz + z * t + # TODO: REMOVABLE? + # representative point to calculate if obstacles are in front or behind + # if point_on_surface is None: + # self.belonging_point = self.polygon.point_on_surface() + # else: + # self.belonging_point = point_on_surface + # self.belonging_point = np.array(self.belonging_point) + # internal 2D coordinates-system to translate projections matrix + # only defined if needed later on + self._projection_matrix = None + self._projected_polygon = None + + def get_3D_shades_from(self, solar_zenith, solar_azimuth, *others): + # project ``others`` onto this shaded surface + # return the shade shapely object + solar_vec = solar_vector(solar_zenith, solar_azimuth) # x,y,z + normal_vec = normal_vector_of_surface(self.tilt, self.azimuth) # a,b,c + + def point_projection(vertex): # vertex := Px, Py, Pz + t = -(normal_vec @ vertex) / (solar_vec @ normal_vec) # Eq. (20) + p_prime = vertex + (t * solar_vec.T).T # Eq. (19) return p_prime - corns = np.array([projection(corner) for corner in self.polygon.coords[:-1]]) + # undo surface rotations to make the third coordinate zero - # TODO - new_points = from_3d_plane_to_2d( - corners_projected, - self.base.surface_tilt, - self.base.surface_azimuth, + _projection_matrix = self._get_projection_matrix_matrix() + _projected_polygon = self._get_self_projected_polygon() + + projected_vertices = np.fromiter( + map(point_projection, other.polygon.coords[:-1]), dtype=(float, 3) ) - # create a polygon - polygon = sp.geometry.Polygon(new_points) - # create a convex hull - hull = ConvexHull(new_points) - # create a shapely polygon - shapely_hull = sp.geometry.Polygon(hull.points[hull.vertices]) - return polygon, shapely_hull + + def get_3D_shade_from_flat_surface(other): + # Section 4.3 in [2] + if projected_vertices.ndim == 1: + projected_vertices = projected_vertices.reshape(1, -1) + vertices_2d = np.fromiter( + map(_projection_matrix.__matmul__, projected_vertices), + dtype=(float, 3), + ) + if not np.allclose(vertices_2d[:, 2], 0.0, atol=1e-10): + raise RuntimeError( + "The third coordinate should be zero, check input " + + "parameters are consistent with the projection plane." + + " If you see this, the error is on me. I fkd up." + ) + vertices_2d = np.delete(vertices_2d, 2, axis=1) + # create a 2D polygon + polygon = sp.Polygon(vertices_2d).intersection(_projected_polygon) + return polygon + + return tuple(map(get_shade_from_flat_surface, others)) + + def _get_projection_matrix_matrix(self): + if self._projection_matrix is None: + self._projection_matrix = (Rz(90 + self.azimuth) @ Rx(self.tilt)).T + return self._projection_matrix + + def _get_self_projected_polygon(self): + if self._projected_polygon is None: + _projection_matrix = self._get_projection_matrix_matrix() + self._projected_polygon = sp.Polygon( + ( + _projection_matrix @ vertex + for vertex in self.polygon.coords[:-1] + ) + ) + return self._projected_polygon -class RectangularSurface(CoplanarSurface): +# %% +class RectangularSurface(FlatSurface): def __init__( self, center, surface_azimuth, surface_tilt, axis_tilt, width, length ): @@ -328,6 +374,8 @@ def __init__( width: width of the surface length: length of the surface TODO: which one is perpendicular to azimuth? + + """ self.center = np.array(center) x_c, y_c, z_c = center @@ -345,25 +393,103 @@ def __init__( ] ) # rotate corners to match the surface orientation - rX, rY, rZ = Rx(row_tilt), Ry(row_azimuth), Rz(axis_tilt) - rot = rX @ rY @ rZ - self.corners = np.array( - [rot @ (corner - self.center) + self.center for corner in corners] - ) + # note pvlib convention uses a left-handed azimuth rotation + rot = Rz(180 - surface_azimuth) @ Ry(axis_tilt) @ Rx(surface_tilt) self.shapely_obj = sp.LinearRing( [rot @ (corner - self.center) + self.center for corner in corners] ) + tilt, azimuth = self._calc_surface_tilt_and_azimuth(rot) + super().__init__(azimuth, tilt, self.shapely_obj) + self._projection_matrix = rot.T + + @classmethod + def _calc_surface_tilt_and_azimuth(cls, rotation_matrix): + # tz as in K. Anderson and M. Mikofski paper, somefig + tz_x, tz_y, tz_z = rotation_matrix[:, 2] # := rot @ [0, 0, 1].T + tilt = acosd(tz_z) + azimuth = atan2d(tz_y, tz_x) + return tilt, azimuth + + def plot(self, ax=None, **kwargs): + if ax is None: + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + x, y, z = np.hsplit( + np.array(self.shapely_obj.coords[:-1]).flatten(order="F"), 3 + ) + ax.plot_trisurf(x, y, z, triangles=((0, 1, 2), (0, 2, 3)), **kwargs) + return ax - # @classmethod - # def from_slope_values(cls, center, slope_azimuth, slope_tilt, width, length): - # return cls(center, slope_azimuth, slope_tilt, width, length) - # @classmethod - # def from_Array_values(cls, center, surface_azimuth, surface_tilt, axis_tilt, width, length): - # return cls(center, surface_azimuth, surface_tilt, axis_tilt, width, length) +# %% +# 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 | [°] | +# +------------------------+----------+-------------+----------+-------+ +# +# + +# 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) +) +solar_azimuth = 120 +solar_zenith = [60, 80] + +# %% +# Fixed Tilt +# ---------- + +fig = plt.figure() +ax = fig.add_subplot(111, projection="3d") +field = RectangularSurface([20 / 2, 10 / 2, 0], 0, 0, 0, 10, 20) +field.plot(ax=ax, color="forestgreen", alpha=0.5) +pv_row1 = RectangularSurface([20 / 2, 0, 2 / 2], 180, 90, 0, 2, 20) +pv_row1.plot(ax=ax, color="darkblue", alpha=0.5) +pv_row2 = RectangularSurface([20 / 2, 10, 2 / 2], 180, 90, 0, 2, 20) +pv_row2.plot(ax=ax, color="darkblue", alpha=0.5) - +# %% +shades = field.get_shades_by(solar_zenith, solar_azimuth, pv_row1, pv_row2) # %% From 732e1f8000a722cf3ee179a8cdf9e11a818a0e33 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sun, 16 Jun 2024 02:27:20 +0200 Subject: [PATCH 08/25] Delete PAR_direct_and_diffuse_sfs.py --- PAR_direct_and_diffuse_sfs.py | 495 ---------------------------------- 1 file changed, 495 deletions(-) delete mode 100644 PAR_direct_and_diffuse_sfs.py diff --git a/PAR_direct_and_diffuse_sfs.py b/PAR_direct_and_diffuse_sfs.py deleted file mode 100644 index e55b6d1b21..0000000000 --- a/PAR_direct_and_diffuse_sfs.py +++ /dev/null @@ -1,495 +0,0 @@ -# %% -import matplotlib.pyplot as plt -import numpy as np -from scipy.spatial import ConvexHull -from scipy.spatial.transform import Rotation as R -from pvlib.tools import sind, cosd - -from itertools import repeat - -# %% -x0, y0, z0 = 0, 0, 1 / 2 # m -W, L = 1, 2 # m -corners0 = np.array( - [ - [x0, y0, z0], - [x0, y0 + W, z0], - [x0 + L, y0, z0], - [x0 + L, y0 + W, z0], - ] -) - - -# %% -# def rotation -def Rz(theta): - return np.array( - [ - [cosd(theta), -sind(theta), 0], - [sind(theta), cosd(theta), 0], - [0, 0, 1], - ] - ) - - -def Ry(theta): - return np.array( - [ - [cosd(theta), 0, sind(theta)], - [0, 1, 0], - [-sind(theta), 0, cosd(theta)], - ] - ) - - -def Rx(theta): - return np.array( - [ - [1, 0, 0], - [0, cosd(theta), -sind(theta)], - [0, sind(theta), cosd(theta)], - ] - ) - - -# %% -# rotate corners, general abstraction of row -def rotate_row_with_center_pivot(corners, theta_x, theta_y): - Rot_x, Rot_y = Rx(theta_x), Ry(theta_y) - midpoint = np.mean(corners, axis=0) - Rot = Rot_x @ Rot_y - return np.array( - [Rot @ (corner - midpoint) + midpoint for corner in corners] - ) - - -# %% -# unittest previous function -def test_rotate_row_with_center_pivot(): - corners = np.array( - [ - [0, 0, 0], - [0, 1, 0], - [1, 0, 0], - [1, 1, 0], - ] - ) - corners_rotated = rotate_row_with_center_pivot(corners, -90, 0) - print(corners_rotated) - - -test_rotate_row_with_center_pivot() - - -# %% -def plot_corners(corners, fig=None, ax=None, **kwargs) -> plt.Figure: - if fig is None: - fig = plt.figure() - if ax is None: - ax = fig.add_subplot(111, projection="3d") - if not isinstance(corners, np.ndarray): - corners = np.array(corners) - x_ = corners[:, 0].reshape(2, -1) - y_ = corners[:, 1].reshape(2, -1) - z_ = corners[:, 2].reshape(2, -1) - ax.plot_surface(x_, y_, z_, **kwargs) - return fig - - -# %% -fig = plot_corners(corners0, alpha=0.5) -fig.show() - -# %% -rotated_fixed_tilt = rotate_row_with_center_pivot(corners0, -90, 0) -fig = plot_corners(rotated_fixed_tilt, alpha=0.5) -fig.show() - - -# %% -def solar_vector(zenith, azimuth): - # Eq. (8), but in terms of zenith instead of elevation - return np.array( - [ - sind(zenith) * sind(azimuth), - sind(zenith) * cosd(azimuth), - cosd(zenith), - ] - ) - - -def normal_vector_of_row(tilt, azimuth): - # Eq. (18) - return np.array( - [ - sind(tilt) * cosd(azimuth), - sind(tilt) * sind(azimuth), - cosd(tilt), - ] - ) - - -def corners_projection_onto_plane( - corners, solar_zenith, solar_azimuth, plane_tilt, plane_azimuth -): - x, y, z = solar_vector(solar_zenith, solar_azimuth) - a, b, c = normal_vector_of_row(plane_tilt, plane_azimuth) - - def projection(corner): - Px, Py, Pz = corner - # Eq. (20) - t = -(a * Px + b * Py + c * Pz) / (a * x + b * y + c * z) - # Eq. (19) - p_prime = np.array([Px + x * t, Py + y * t, Pz + z * t]) - return p_prime - - return np.array([projection(corner) for corner in corners]) - - -# %% -# set hypothetical values -solar_zenith = 30 -solar_azimuth = 180 -row_tilt = 30 -row_azimuth = 180 -plane_tilt = -10 -plane_azimuth = 90 - -# create the row -corners1 = rotate_row_with_center_pivot(corners0, row_tilt, 0) -corners_projected = corners_projection_onto_plane( - corners1, solar_zenith, solar_azimuth, plane_tilt, plane_azimuth -) - -fig = plt.figure() -ax = fig.add_subplot(111, projection="3d") -ax.set_xlim(0, 2) -ax.set_ylim(0, 3) -ax.set_zlim(0, 1) -_ = plot_corners(corners1, ax=ax, alpha=0.5, color="grey") -_ = plot_corners(corners_projected, ax=ax, alpha=0.5, color="brown") -fig.show() - -# %% -# Compare corners_projected with corners after rotation to new coordinate systems of the projection plane, to remove the third coordinate - - -def from_3d_plane_to_2d(points, tilt, azimuth): - # Section 4.3 in [2] - R_z = Rz(90 + azimuth) - R_x = Rx(tilt) - rot = (R_z @ R_x).T - if points.ndim == 1: - points = points.reshape(1, -1) - new_points = np.fromiter(map(rot.__matmul__, points), dtype=(float, 3)) - assert np.allclose( - new_points[:, 2], 0 - ), "The third coordinate should be zero, check input parameters are consistent with the projection plane." - return np.delete(new_points, 2, axis=1) - - -new_points = from_3d_plane_to_2d(corners_projected, plane_tilt, plane_azimuth) - -print(new_points) - -# %% -# -########################## -# OOP Implementation -########################## -import numpy as np -import pandas as pd -import shapely as sp -from pvlib.tools import sind, cosd, acosd -import matplotlib.pyplot as plt - -from abc import ABC, abstractmethod - - -def atan2d(y, x): - return np.degrees(np.arctan2(y, x)) - - -def Rz(theta): - return np.array( - [ - [cosd(theta), -sind(theta), 0], - [sind(theta), cosd(theta), 0], - [0, 0, 1], - ] - ) - - -def Ry(theta): - return np.array( - [ - [cosd(theta), 0, sind(theta)], - [0, 1, 0], - [-sind(theta), 0, cosd(theta)], - ] - ) - - -def Rx(theta): - return np.array( - [ - [1, 0, 0], - [0, cosd(theta), -sind(theta)], - [0, sind(theta), cosd(theta)], - ] - ) - - -def solar_vector(zenith, azimuth): - # Eq. (8), but in terms of zenith instead of elevation - return np.array( - [ - sind(zenith) * sind(azimuth), - sind(zenith) * cosd(azimuth), - cosd(zenith), - ] - ).T - - -def normal_vector_of_surface(tilt, azimuth): - # Eq. (18) - return np.array( - [ - sind(tilt) * cosd(azimuth), - sind(tilt) * sind(azimuth), - cosd(tilt), - ] - ).T - - -# %% -class FlatSurface: - def __init__( - self, azimuth, tilt, polygon_boundaries - ): - """ - - .. warning:: - - This constructor does not make any check on the input parameters. - Use any of the class methods ``from_*`` to create an object. - - Parameters - ---------- - azimuth : float - Surface azimuth, angle at which it points downwards. 0° is North, 90° is East, 180° is South, 270° is West. In degrees [°]. - tilt : float - Surface tilt, angle at which it is inclined with respect to the horizontal plane. Tilted downwards ``azimuth``. 0° is horizontal, 90° is vertical. In degrees [°]. - polygon : shapely.LinearRing or array[N, 3] - Shapely Polygon or boundaries to build it. If you need to specify holes, provide your custom LinearRing. - """ - self.azimuth = azimuth - self.tilt = tilt - # works for polygon_boundaries := array[N, 3] | shapely.Polygon - self.polygon = sp.LinearRing(polygon_boundaries) - # TODO: REMOVABLE? - # representative point to calculate if obstacles are in front or behind - # if point_on_surface is None: - # self.belonging_point = self.polygon.point_on_surface() - # else: - # self.belonging_point = point_on_surface - # self.belonging_point = np.array(self.belonging_point) - # internal 2D coordinates-system to translate projections matrix - # only defined if needed later on - self._projection_matrix = None - self._projected_polygon = None - - def get_3D_shades_from(self, solar_zenith, solar_azimuth, *others): - # project ``others`` onto this shaded surface - # return the shade shapely object - solar_vec = solar_vector(solar_zenith, solar_azimuth) # x,y,z - normal_vec = normal_vector_of_surface(self.tilt, self.azimuth) # a,b,c - - def point_projection(vertex): # vertex := Px, Py, Pz - t = -(normal_vec @ vertex) / (solar_vec @ normal_vec) # Eq. (20) - p_prime = vertex + (t * solar_vec.T).T # Eq. (19) - return p_prime - - # undo surface rotations to make the third coordinate zero - _projection_matrix = self._get_projection_matrix_matrix() - _projected_polygon = self._get_self_projected_polygon() - - projected_vertices = np.fromiter( - map(point_projection, other.polygon.coords[:-1]), dtype=(float, 3) - ) - - def get_3D_shade_from_flat_surface(other): - # Section 4.3 in [2] - if projected_vertices.ndim == 1: - projected_vertices = projected_vertices.reshape(1, -1) - vertices_2d = np.fromiter( - map(_projection_matrix.__matmul__, projected_vertices), - dtype=(float, 3), - ) - if not np.allclose(vertices_2d[:, 2], 0.0, atol=1e-10): - raise RuntimeError( - "The third coordinate should be zero, check input " - + "parameters are consistent with the projection plane." - + " If you see this, the error is on me. I fkd up." - ) - vertices_2d = np.delete(vertices_2d, 2, axis=1) - # create a 2D polygon - polygon = sp.Polygon(vertices_2d).intersection(_projected_polygon) - return polygon - - return tuple(map(get_shade_from_flat_surface, others)) - - def _get_projection_matrix_matrix(self): - if self._projection_matrix is None: - self._projection_matrix = (Rz(90 + self.azimuth) @ Rx(self.tilt)).T - return self._projection_matrix - - def _get_self_projected_polygon(self): - if self._projected_polygon is None: - _projection_matrix = self._get_projection_matrix_matrix() - self._projected_polygon = sp.Polygon( - ( - _projection_matrix @ vertex - for vertex in self.polygon.coords[:-1] - ) - ) - return self._projected_polygon - - -# %% -class RectangularSurface(FlatSurface): - def __init__( - self, center, surface_azimuth, surface_tilt, axis_tilt, width, length - ): - """ - - .. warning:: - - This constructor does not make any check on the input parameters. - Use any of the class methods ``from_*`` to create an object. - - center: center of the surface - surface_azimuth: azimuth of the surface - surface_tilt: tilt of the surface - width: width of the surface - length: length of the surface - TODO: which one is perpendicular to azimuth? - - - """ - self.center = np.array(center) - x_c, y_c, z_c = center - self.surface_azimuth = surface_azimuth - self.surface_tilt = surface_tilt - self.axis_tilt = axis_tilt - self.width = width - self.length = length - corners = np.array( - [ - [x_c - length / 2, y_c - width / 2, z_c], - [x_c - length / 2, y_c + width / 2, z_c], - [x_c + length / 2, y_c + width / 2, z_c], - [x_c + length / 2, y_c - width / 2, z_c], - ] - ) - # rotate corners to match the surface orientation - # note pvlib convention uses a left-handed azimuth rotation - rot = Rz(180 - surface_azimuth) @ Ry(axis_tilt) @ Rx(surface_tilt) - self.shapely_obj = sp.LinearRing( - [rot @ (corner - self.center) + self.center for corner in corners] - ) - tilt, azimuth = self._calc_surface_tilt_and_azimuth(rot) - super().__init__(azimuth, tilt, self.shapely_obj) - self._projection_matrix = rot.T - - @classmethod - def _calc_surface_tilt_and_azimuth(cls, rotation_matrix): - # tz as in K. Anderson and M. Mikofski paper, somefig - tz_x, tz_y, tz_z = rotation_matrix[:, 2] # := rot @ [0, 0, 1].T - tilt = acosd(tz_z) - azimuth = atan2d(tz_y, tz_x) - return tilt, azimuth - - def plot(self, ax=None, **kwargs): - if ax is None: - fig = plt.figure() - ax = fig.add_subplot(111, projection="3d") - x, y, z = np.hsplit( - np.array(self.shapely_obj.coords[:-1]).flatten(order="F"), 3 - ) - ax.plot_trisurf(x, y, z, triangles=((0, 1, 2), (0, 2, 3)), **kwargs) - return ax - - -# %% -# 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 | [°] | -# +------------------------+----------+-------------+----------+-------+ -# -# - -# 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) -) - -solar_azimuth = 120 -solar_zenith = [60, 80] - -# %% -# Fixed Tilt -# ---------- - -fig = plt.figure() -ax = fig.add_subplot(111, projection="3d") -field = RectangularSurface([20 / 2, 10 / 2, 0], 0, 0, 0, 10, 20) -field.plot(ax=ax, color="forestgreen", alpha=0.5) -pv_row1 = RectangularSurface([20 / 2, 0, 2 / 2], 180, 90, 0, 2, 20) -pv_row1.plot(ax=ax, color="darkblue", alpha=0.5) -pv_row2 = RectangularSurface([20 / 2, 10, 2 / 2], 180, 90, 0, 2, 20) -pv_row2.plot(ax=ax, color="darkblue", alpha=0.5) - - - -# %% -shades = field.get_shades_by(solar_zenith, solar_azimuth, pv_row1, pv_row2) -# %% From 53d03b52cab2886a1554a91aaaf5537ebca0833d Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Tue, 18 Jun 2024 03:26:18 +0200 Subject: [PATCH 09/25] Spatial module and examples beginning --- docs/examples/agrivoltaics/README.rst | 2 + .../plot_par_direct_shading0_fixed_tilt.py | 176 ++++++++ .../plot_spatial_row_to_row_shading.py | 75 ++++ pvlib/__init__.py | 1 + pvlib/spatial.py | 411 ++++++++++++++++++ pvlib/tests/test_spatial.py | 178 ++++++++ 6 files changed, 843 insertions(+) create mode 100644 docs/examples/agrivoltaics/README.rst create mode 100644 docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py create mode 100644 docs/examples/shading/plot_spatial_row_to_row_shading.py create mode 100644 pvlib/spatial.py create mode 100644 pvlib/tests/test_spatial.py 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..4e2fdbdcb6 --- /dev/null +++ b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py @@ -0,0 +1,176 @@ +""" +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 shapely.plotting +import pandas as pd +import numpy as np +import shapely as sp +from pvlib.spatial import RectangularSurface + +# 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) +) +solar_azimuth = 155 +solar_zenith = 70 + +# %% +# Fixed Tilt +# ---------- + +field = RectangularSurface( # crops top surface + center=[10, 5, 0], + azimuth=0, + tilt=0, + axis_tilt=0, + width=10, + length=20, +) +pv_row1 = RectangularSurface( # south-most row (lowest Y-coordinate) + center=[10, -10 / 2 + 5, 2 / 2], + azimuth=180, + tilt=90, + axis_tilt=0, + width=2, + length=20, +) +pv_row2 = RectangularSurface( # north-most row (highest Y-coordinate) + center=[10, 10 / 2 + 5, 2 / 2], + azimuth=180, + tilt=90, + axis_tilt=0, + width=2, + length=20, +) + +# %% +shades_3d = field.get_3D_shades_from( + solar_zenith, solar_azimuth, pv_row1, pv_row2 +) +print(shades_3d) +# %% +shades_2d = field.get_2D_shades_from( + solar_zenith, solar_azimuth, pv_row1, pv_row2 +) +print(shades_2d) + +# %% +# Plot both the 3D and 2D shades +# ------------------------------ + +field_style = {"color": "forestgreen", "alpha": 0.5} +row_style = {"color": "darkblue", "alpha": 0.5} +shade_style = {"color": "dimgrey", "alpha": 0.8} + +fig = plt.figure(figsize=(10, 10)) +gs = fig.add_gridspec(10, 1) + +ax1 = fig.add_subplot(gs[0:6, 0], projection="3d") +ax2 = fig.add_subplot(gs[8:, 0]) + +ax1.view_init(elev=45, azim=-45) +field.plot(ax=ax1, **field_style) +pv_row1.plot(ax=ax1, **row_style) +pv_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]") + +field2d = field.representation_in_2D_space() +field_style_2d = {**field_style, "add_points": False} +shapely.plotting.plot_polygon(field2d, ax=ax2, **field_style_2d) +shade_style_2d = {**shade_style, "add_points": False} +for shade in shades_2d.geoms: + shapely.plotting.plot_polygon(shade, ax=ax2, **shade_style_2d) + +# %% +beam_shaded_fraction = ( + sum(shade.area for shade in shades_2d.geoms) / field2d.area +) +print(beam_shaded_fraction) + +# %% +# 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..956771dd5d --- /dev/null +++ b/docs/examples/shading/plot_spatial_row_to_row_shading.py @@ -0,0 +1,75 @@ +""" +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 shading fraction in 3D and 2D space with the help of the +# :py:mod:`~pvlib.spatial` module and its classes. +# +# First section of the example will show how to calculate the shading fraction +# for an instantaneous time. The second section will show how to calculate the +# shading fraction for a time range. + +from pvlib.spatial import RectangularSurface +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d.art3d import Poly3DCollection +import shapely.plotting + +# %% +# 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. + +solar_azimuth = 180 # degrees +solar_zenith = 70 # degrees + +# Create the first row of panels +row1 = RectangularSurface( # south-most row + center=[0, 0, 3], azimuth=180, tilt=20, axis_tilt=0, width=2, length=20 +) + +# Create the second row of panels +row2 = RectangularSurface( # north-most row + center=[0, 3, 3], azimuth=180, tilt=50, axis_tilt=0, width=2, length=20 +) + +# Calculate the shadow +shades_3d = row2.get_3D_shades_from(solar_zenith, solar_azimuth, row1) +shades_2d = row2.get_2D_shades_from(solar_zenith, solar_azimuth, row1) + +# %% +# Plot the scene and the shadow +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +row_style = {"color": "darkblue", "alpha": 0.5} +shade_style = {"color": "dimgrey", "alpha": 0.8} + +fig = plt.figure(figsize=(10, 10)) +gs = fig.add_gridspec(10, 1) + +ax1 = fig.add_subplot(gs[0:6, 0], projection="3d") +ax2 = fig.add_subplot(gs[8:, 0]) + +ax1.view_init(elev=45, azim=-45) +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)) + +row_style_2d = {**row_style, "add_points": False} +row2_2d = row2.representation_in_2D_space(solar_zenith, solar_azimuth) +print(f"{row2_2d=}") +shapely.plotting.plot_polygon(row2_2d, ax=ax2, **row_style_2d) +shade_style_2d = {**shade_style, "add_points": False} +for shade in shades_2d.geoms: + shapely.plotting.plot_polygon(shade, ax=ax2, **shade_style_2d) 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..06e9973f1f --- /dev/null +++ b/pvlib/spatial.py @@ -0,0 +1,411 @@ +""" +Spatial functions for shading and 3D scenes analysis. +""" + +from pvlib.tools import sind, cosd, acosd, atan2d +import numpy as np +import shapely as sp +import matplotlib.pyplot as plt +from scipy.spatial.transform import Rotation + + +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. + + + 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 polygon(self): + return self._polygon + + def __init__(self, azimuth, tilt, polygon_boundaries): + # Wrap these two attributes in 2D ndarrays for consistency + self._azimuth = np.array(azimuth) + self._tilt = np.array(tilt) + # works for polygon_boundaries := array[N, 3] | shapely.Polygon + self._polygon = sp.Polygon(polygon_boundaries) + # internal 2D coordinates-system to translate projections matrix + # only defined if needed later on + self._rotation_inverse = None + self._projected_polygon = None + + 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 + ------- + tuple[shapely.Polygon] + Shapely Polygon objects representing the shades on this surface. + """ + 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 + ) + + plane_point = np.array( + self._polygon.exterior.coords[0] + ) # any point on the plane + + 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 = ((plane_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 + shade = sp.Polygon(projected_vertices).intersection( + _polygon, grid_size=1e-12 + ) + return shade if isinstance(shade, sp.Polygon) else sp.Polygon() + + return sp.MultiPolygon(map(get_3D_shade_from_flat_surface, others)) + + def get_2D_shades_from(self, solar_zenith, solar_azimuth, *others): + 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_origin_plane(vertex): # vertex -> Px, Py, Pz + # Eq. (20), projects to plane that goes through the origin (0,0,0) + t = -(vertex @ normal_vec) / (solar_vec @ normal_vec) + p_prime = vertex + (t * solar_vec.T).T # Eq. (19) + return p_prime + + # undo surface rotations to make the third coordinate zero + _projection = self._get_projection_transform() + # and clip the 2D shades to the surface in 2D + # _self_projected_polygon = self.representation_in_2D_space() + # print(f"{_self_projected_polygon=}") + + # Section 4.3 in [2] + def transform_to_2D_reference_plane(vertices): + vertices_2d = _projection.apply(vertices) + if not np.allclose(vertices_2d[:, 2], 0.0, atol=1e-10): + raise RuntimeError( + "Non-null third coordinate in 2D projection!" + ) # for debugging purposes; TODO: remove <<<<<< !!!!!! ###### + return np.delete(vertices_2d, 2, axis=1) + + def get_2D_shade_from_flat_surface(other): + coords_to_project = np.array(other.polygon.exterior.coords[:-1]) + projected_vertices = np.fromiter( + map(project_point_to_origin_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 + projected_vertices_2d = transform_to_2D_reference_plane( + projected_vertices + ) + shade = sp.Polygon(projected_vertices_2d) + # .intersection( + # _self_projected_polygon + # ) + print(f"{sp.Polygon(projected_vertices_2d)=}") + print(f"{shade=}") + return shade if isinstance(shade, sp.Polygon) else sp.Polygon() + + return sp.MultiPolygon(map(get_2D_shade_from_flat_surface, others)) + + def _get_projection_transform(self): + if self._rotation_inverse is None: + self._rotation_inverse = Rotation.from_euler( + "ZXZ", [-self._azimuth, -self._tilt, 0], degrees=True + ).inv() + return self._rotation_inverse + + def representation_in_2D_space(self, solar_zenith, solar_azimuth): + if self._projected_polygon is None: + _projection = self._get_projection_transform() + 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_origin_plane(vertex): # vertex -> Px, Py, Pz + # Eq. (20), projects to plane that goes through the origin (0,0,0) + t = -(vertex @ normal_vec) / (solar_vec @ normal_vec) + p_prime = vertex + (t * solar_vec.T).T # Eq. (19) + return p_prime + vertices_to_project = self._polygon.exterior.coords[:-1] + origin_plane_vertices = np.fromiter( + map(project_point_to_origin_plane, vertices_to_project), + dtype=(float, 3), + count=len(vertices_to_project), # speeds up allocation + ) + projected_vertices = _projection.apply(origin_plane_vertices) + self._projected_polygon = sp.Polygon( + np.delete(projected_vertices, 2, axis=1) + ) + 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. + """ + return sp.ops.unary_union(shades) + + def plot(self, ax=None, **kwargs): + pass # TODO: implement this method + + +# %% +class RectangularSurface(FlatSurface): + def __init__(self, center, azimuth, tilt, axis_tilt, width, length): + """ + 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:`FlatSurface` handy for rectangular surfaces like PV arrays. + + See :py:class:`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 [°]. + 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 + """ + self.center = np.array(center) + 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 + _rotation = Rotation.from_euler( + "ZXZ", [-azimuth, -tilt, axis_tilt], degrees=True + ) + self._polygon = sp.Polygon(_rotation.apply(corners) + center) + super().__init__(azimuth, tilt, self._polygon) + self._rotation_inverse = _rotation.inv() + + @classmethod + def _calc_surface_tilt_and_azimuth(cls, rotation_matrix): + """ + Given the rotation matrix that results from the surface orientation in + terms of an ``surface_azimuth``, ``surface_tilt`` and ``axis_tilt``, + calculate the resulting tilt and azimuth angles of the surface. + + Parameters + ---------- + rotation_matrix : array[3, 3] + Rotation matrix. + + Returns + ------- + tilt, azimuth : float, float + Surface tilt and azimuth angles in degrees. + """ + # tz as in K. Anderson and M. Mikofski paper, Fig. 1 + tz_x, tz_y, tz_z = rotation_matrix[:, 2] # := rot @ [0, 0, 1].T + tilt = acosd(tz_z) + azimuth = atan2d(tz_y, tz_x) + return tilt, azimuth + + def plot(self, ax=None, **kwargs): + """ + Plot the rectangular surface. + + Parameters + ---------- + ax : matplotlib.axes.Axes, 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 + ) + ax.plot_trisurf(x, y, z, triangles=((0, 1, 2), (0, 2, 3)), **kwargs) + return ax diff --git a/pvlib/tests/test_spatial.py b/pvlib/tests/test_spatial.py new file mode 100644 index 0000000000..312c741717 --- /dev/null +++ b/pvlib/tests/test_spatial.py @@ -0,0 +1,178 @@ +"""Tests spatial submodule.""" + +from pvlib import spatial +from numpy.testing import assert_allclose +import shapely + +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 = shapely.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, shapely.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 == shapely.Polygon(polygon) + assert isinstance(surface.polygon, shapely.Polygon) + + +def test_FlatSurface_readonly_properties(): + surface_azimuth = 180 + surface_tilt = 0 + polygon = shapely.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.center == center + assert surface.azimuth == surface_azimuth + assert surface.tilt == surface_tilt + assert surface.axis_tilt == axis_tilt + assert surface.width == width + assert surface.length == length + # 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.center == center + assert surface.azimuth == surface_azimuth + assert surface.tilt == surface_tilt + assert surface.axis_tilt == axis_tilt + assert surface.width == width + assert surface.length == length + + +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 \ No newline at end of file From d83d6aa9691b672eee13ce90a9e9f0cbeb21149f Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Thu, 20 Jun 2024 02:24:32 +0200 Subject: [PATCH 10/25] Fix rotations, refactor --- .../plot_par_direct_shading0_fixed_tilt.py | 8 +- .../plot_spatial_row_to_row_shading.py | 13 +- pvlib/spatial.py | 222 +++++++++--------- 3 files changed, 122 insertions(+), 121 deletions(-) diff --git a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py index 4e2fdbdcb6..0601d20f2f 100644 --- a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py +++ b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py @@ -54,8 +54,6 @@ from mpl_toolkits.mplot3d.art3d import Poly3DCollection import shapely.plotting import pandas as pd -import numpy as np -import shapely as sp from pvlib.spatial import RectangularSurface # Kärrbo Prästgård, Västerås, Sweden @@ -79,7 +77,7 @@ field = RectangularSurface( # crops top surface center=[10, 5, 0], - azimuth=0, + azimuth=180, # chosen instead of 0 (north) for intuitive visualization tilt=0, axis_tilt=0, width=10, @@ -141,8 +139,8 @@ ax1.axis("equal") ax1.set_zlim(0) -ax1.set_xlabel("West(-) <> East(+) [m]") -ax1.set_ylabel("South(-) <> North(+) [m]") +ax1.set_xlabel("West(-) East(+) [m]") +ax1.set_ylabel("South(-) North(+) [m]") field2d = field.representation_in_2D_space() field_style_2d = {**field_style, "add_points": False} diff --git a/docs/examples/shading/plot_spatial_row_to_row_shading.py b/docs/examples/shading/plot_spatial_row_to_row_shading.py index 956771dd5d..470b58c819 100644 --- a/docs/examples/shading/plot_spatial_row_to_row_shading.py +++ b/docs/examples/shading/plot_spatial_row_to_row_shading.py @@ -26,16 +26,16 @@ # will be 30 degrees. The distance between the two rows will be 3 meters. solar_azimuth = 180 # degrees -solar_zenith = 70 # degrees +solar_zenith = 75 # degrees # Create the first row of panels row1 = RectangularSurface( # south-most row - center=[0, 0, 3], azimuth=180, tilt=20, axis_tilt=0, width=2, length=20 + center=[0, 0, 3], azimuth=165, tilt=20, axis_tilt=10, width=2, length=20 ) # Create the second row of panels row2 = RectangularSurface( # north-most row - center=[0, 3, 3], azimuth=180, tilt=50, axis_tilt=0, width=2, length=20 + center=[0, 3, 3], azimuth=165, tilt=20, axis_tilt=10, width=2, length=20 ) # Calculate the shadow @@ -66,8 +66,13 @@ 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]") + row_style_2d = {**row_style, "add_points": False} -row2_2d = row2.representation_in_2D_space(solar_zenith, solar_azimuth) +row2_2d = row2.representation_in_2D_space() print(f"{row2_2d=}") shapely.plotting.plot_polygon(row2_2d, ax=ax2, **row_style_2d) shade_style_2d = {**shade_style, "add_points": False} diff --git a/pvlib/spatial.py b/pvlib/spatial.py index 06e9973f1f..6e550996ba 100644 --- a/pvlib/spatial.py +++ b/pvlib/spatial.py @@ -2,7 +2,7 @@ Spatial functions for shading and 3D scenes analysis. """ -from pvlib.tools import sind, cosd, acosd, atan2d +from pvlib.tools import sind, cosd import numpy as np import shapely as sp import matplotlib.pyplot as plt @@ -33,9 +33,10 @@ def _solar_vector(zenith, azimuth): .. [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`. + .. [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( @@ -75,9 +76,10 @@ def _plane_normal_vector(tilt, azimuth): 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`. + .. [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( @@ -85,7 +87,6 @@ def _plane_normal_vector(tilt, azimuth): ) -# %% class FlatSurface: """ Represents a flat surface in 3D space with a given azimuth and tilt and @@ -113,7 +114,12 @@ class FlatSurface: polygon : shapely.Polygon or array[N, 3] Shapely Polygon or boundaries to build it. Holes are ignored for now. - + normal_rotation : 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 ---------- @@ -126,9 +132,9 @@ class FlatSurface: 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`. + backtracking in single-axis trackers on rolling terrain. J. Renewable + Sustainable Energy 1 March 2024; 16 (2): 023504. + :doi:`10.1063/5.0202220`. """ @property @@ -139,21 +145,70 @@ def azimuth(self): def tilt(self): return self._tilt + @property + def normal_rotation(self): + return self._normal_rotation + @property def polygon(self): return self._polygon - def __init__(self, azimuth, tilt, polygon_boundaries): - # Wrap these two attributes in 2D ndarrays for consistency + @property + def reference_point(self): + return self._reference_point + + def __init__( + self, + azimuth, + tilt, + polygon_boundaries, + normal_rotation=0.0, + reference_point=None, + ): self._azimuth = np.array(azimuth) self._tilt = np.array(tilt) + self._normal_rotation = np.array(normal_rotation) # works for polygon_boundaries := array[N, 3] | shapely.Polygon self._polygon = sp.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. + + Parameters + ---------- + point : array-like, shape (N, 3) + Point to project. + + Returns + ------- + array-like, shape (N, 2) + Projected point. + """ + # undo surface rotations to make the third coordinate zero + _projection = self._get_projection_transform() + # Section 4.3 in [2] + 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. @@ -176,8 +231,9 @@ def get_3D_shades_from(self, solar_zenith, solar_azimuth, *others): Returns ------- - tuple[shapely.Polygon] + shapely.MultiPolygon Shapely Polygon objects representing the shades on this surface. + These shades may overlap. """ solar_vec = solar_vec = _solar_vector( # Eq. (8) -> x,y,z solar_zenith, solar_azimuth @@ -186,14 +242,10 @@ def get_3D_shades_from(self, solar_zenith, solar_azimuth, *others): self._tilt, self._azimuth ) - plane_point = np.array( - self._polygon.exterior.coords[0] - ) # any point on the plane - 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 = ((plane_point - vertex) @ normal_vec) / ( + t = ((self._reference_point - vertex) @ normal_vec) / ( solar_vec @ normal_vec ) p_prime = vertex + (t * solar_vec.T).T # Eq. (19) @@ -217,86 +269,49 @@ def get_3D_shade_from_flat_surface(other): return sp.MultiPolygon(map(get_3D_shade_from_flat_surface, others)) def get_2D_shades_from(self, solar_zenith, solar_azimuth, *others): - 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 - ) + shades = self.get_3D_shades_from(solar_zenith, solar_azimuth, *others) - def project_point_to_origin_plane(vertex): # vertex -> Px, Py, Pz - # Eq. (20), projects to plane that goes through the origin (0,0,0) - t = -(vertex @ normal_vec) / (solar_vec @ normal_vec) - p_prime = vertex + (t * solar_vec.T).T # Eq. (19) - return p_prime - - # undo surface rotations to make the third coordinate zero - _projection = self._get_projection_transform() # and clip the 2D shades to the surface in 2D - # _self_projected_polygon = self.representation_in_2D_space() - # print(f"{_self_projected_polygon=}") - - # Section 4.3 in [2] - def transform_to_2D_reference_plane(vertices): - vertices_2d = _projection.apply(vertices) - if not np.allclose(vertices_2d[:, 2], 0.0, atol=1e-10): - raise RuntimeError( - "Non-null third coordinate in 2D projection!" - ) # for debugging purposes; TODO: remove <<<<<< !!!!!! ###### - return np.delete(vertices_2d, 2, axis=1) + _self_projected_polygon = self.representation_in_2D_space() def get_2D_shade_from_flat_surface(other): - coords_to_project = np.array(other.polygon.exterior.coords[:-1]) - projected_vertices = np.fromiter( - map(project_point_to_origin_plane, coords_to_project), - dtype=(float, 3), - count=len(coords_to_project), # speeds up allocation - ) + coords_to_project = np.array(other.exterior.coords[:-1]) + # create shapely shade object and bound it to the surface - projected_vertices_2d = transform_to_2D_reference_plane( - projected_vertices + projected_vertices_2d = self._transform_to_2D(coords_to_project) + shade = sp.Polygon(projected_vertices_2d).intersection( + _self_projected_polygon ) - shade = sp.Polygon(projected_vertices_2d) - # .intersection( - # _self_projected_polygon - # ) - print(f"{sp.Polygon(projected_vertices_2d)=}") - print(f"{shade=}") + # Return Polygons only, geometries with empty area are omitted return shade if isinstance(shade, sp.Polygon) else sp.Polygon() - return sp.MultiPolygon(map(get_2D_shade_from_flat_surface, others)) + return sp.MultiPolygon( + map(get_2D_shade_from_flat_surface, shades.geoms) + ) def _get_projection_transform(self): if self._rotation_inverse is None: self._rotation_inverse = Rotation.from_euler( - "ZXZ", [-self._azimuth, -self._tilt, 0], degrees=True - ).inv() + "zxz", + [self._azimuth - 180, -self._tilt, -self._normal_rotation], + degrees=True, + ) return self._rotation_inverse - def representation_in_2D_space(self, solar_zenith, solar_azimuth): - if self._projected_polygon is None: - _projection = self._get_projection_transform() - 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 representation_in_2D_space(self): + """ + Get the 2D representation of the surface in its local plane, with + respect to the reference point. - def project_point_to_origin_plane(vertex): # vertex -> Px, Py, Pz - # Eq. (20), projects to plane that goes through the origin (0,0,0) - t = -(vertex @ normal_vec) / (solar_vec @ normal_vec) - p_prime = vertex + (t * solar_vec.T).T # Eq. (19) - return p_prime - vertices_to_project = self._polygon.exterior.coords[:-1] - origin_plane_vertices = np.fromiter( - map(project_point_to_origin_plane, vertices_to_project), - dtype=(float, 3), - count=len(vertices_to_project), # speeds up allocation - ) - projected_vertices = _projection.apply(origin_plane_vertices) + Returns + ------- + shapely.Polygon + 2D representation of the surface. + """ + if self._projected_polygon is None: + vertices = np.array(self._polygon.exterior.coords[:-1]) self._projected_polygon = sp.Polygon( - np.delete(projected_vertices, 2, axis=1) + self._transform_to_2D(vertices) ) return self._projected_polygon @@ -321,7 +336,6 @@ def plot(self, ax=None, **kwargs): pass # TODO: implement this method -# %% class RectangularSurface(FlatSurface): def __init__(self, center, azimuth, tilt, axis_tilt, width, length): """ @@ -348,7 +362,6 @@ def __init__(self, center, azimuth, tilt, axis_tilt, width, length): length: length of the surface Perpendicular to the surface azimuth """ - self.center = np.array(center) corners = np.array( [ [-length / 2, -width / 2, 0], @@ -359,35 +372,20 @@ def __init__(self, center, azimuth, tilt, axis_tilt, width, length): ) # 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 ) - self._polygon = sp.Polygon(_rotation.apply(corners) + center) - super().__init__(azimuth, tilt, self._polygon) - self._rotation_inverse = _rotation.inv() - - @classmethod - def _calc_surface_tilt_and_azimuth(cls, rotation_matrix): - """ - Given the rotation matrix that results from the surface orientation in - terms of an ``surface_azimuth``, ``surface_tilt`` and ``axis_tilt``, - calculate the resulting tilt and azimuth angles of the surface. - - Parameters - ---------- - rotation_matrix : array[3, 3] - Rotation matrix. - - Returns - ------- - tilt, azimuth : float, float - Surface tilt and azimuth angles in degrees. - """ - # tz as in K. Anderson and M. Mikofski paper, Fig. 1 - tz_x, tz_y, tz_z = rotation_matrix[:, 2] # := rot @ [0, 0, 1].T - tilt = acosd(tz_z) - azimuth = atan2d(tz_y, tz_x) - return tilt, azimuth + polygon = sp.Polygon(_rotation.apply(corners) + center) + super().__init__( + azimuth=azimuth, + tilt=tilt, + polygon_boundaries=polygon, + normal_rotation=axis_tilt, + reference_point=center, + ) def plot(self, ax=None, **kwargs): """ From 6750e74a3714e29665c5cc1619229980313558dc Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Thu, 20 Jun 2024 20:13:58 +0200 Subject: [PATCH 11/25] Specially docs and variable naming --- .../plot_par_direct_shading0_fixed_tilt.py | 14 +- docs/sphinx/source/conf.py | 3 + pvlib/spatial.py | 191 ++++++++++++++---- 3 files changed, 167 insertions(+), 41 deletions(-) diff --git a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py index 0601d20f2f..ae9273cc0a 100644 --- a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py +++ b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py @@ -55,6 +55,7 @@ import shapely.plotting import pandas as pd 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 @@ -68,8 +69,11 @@ .union(fall_equinox) .union(winter_solstice) ) -solar_azimuth = 155 -solar_zenith = 70 +solar_position = solarposition.get_solarposition( + dates, latitude, longitude, altitude +) +solar_azimuth = solar_position["azimuth"].iloc[0] +solar_zenith = solar_position["apparent_zenith"].iloc[0] # %% # Fixed Tilt @@ -100,16 +104,12 @@ length=20, ) -# %% shades_3d = field.get_3D_shades_from( solar_zenith, solar_azimuth, pv_row1, pv_row2 ) -print(shades_3d) -# %% shades_2d = field.get_2D_shades_from( - solar_zenith, solar_azimuth, pv_row1, pv_row2 + solar_zenith, solar_azimuth, shades_3d=shades_3d ) -print(shades_2d) # %% # Plot both the 3D and 2D shades 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/pvlib/spatial.py b/pvlib/spatial.py index 6e550996ba..df65412bbc 100644 --- a/pvlib/spatial.py +++ b/pvlib/spatial.py @@ -114,7 +114,7 @@ class FlatSurface: polygon : shapely.Polygon or array[N, 3] Shapely Polygon or boundaries to build it. Holes are ignored for now. - normal_rotation : float, default 0.0° + 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 @@ -146,8 +146,8 @@ def tilt(self): return self._tilt @property - def normal_rotation(self): - return self._normal_rotation + def roll(self): + return self._roll @property def polygon(self): @@ -162,12 +162,12 @@ def __init__( azimuth, tilt, polygon_boundaries, - normal_rotation=0.0, + roll=0.0, reference_point=None, ): self._azimuth = np.array(azimuth) self._tilt = np.array(tilt) - self._normal_rotation = np.array(normal_rotation) + self._roll = np.array(roll) # works for polygon_boundaries := array[N, 3] | shapely.Polygon self._polygon = sp.Polygon(polygon_boundaries) # internal 2D coordinates-system to translate projections matrix @@ -188,6 +188,9 @@ 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) @@ -197,10 +200,23 @@ def _transform_to_2D(self, points): ------- 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 [2] + # 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( @@ -234,6 +250,20 @@ def get_3D_shades_from(self, solar_zenith, solar_azimuth, *others): 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 @@ -268,8 +298,60 @@ def get_3D_shade_from_flat_surface(other): return sp.MultiPolygon(map(get_3D_shade_from_flat_surface, others)) - def get_2D_shades_from(self, solar_zenith, solar_azimuth, *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 = sp.MultiPolygon(shades, shades_3d) # and clip the 2D shades to the surface in 2D _self_projected_polygon = self.representation_in_2D_space() @@ -290,10 +372,29 @@ def get_2D_shade_from_flat_surface(other): ) 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._normal_rotation], + [self._azimuth - 180, -self._tilt, -self._roll], degrees=True, ) return self._rotation_inverse @@ -307,6 +408,11 @@ def representation_in_2D_space(self): ------- 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]) @@ -337,31 +443,48 @@ def plot(self, ax=None, **kwargs): class RectangularSurface(FlatSurface): - def __init__(self, center, azimuth, tilt, axis_tilt, width, length): - """ - 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:`FlatSurface` handy for rectangular surfaces like PV arrays. + """ + 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:`FlatSurface` for information on methods and properties. + 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 [°]. - 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 - """ + 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], @@ -383,8 +506,8 @@ def __init__(self, center, azimuth, tilt, axis_tilt, width, length): azimuth=azimuth, tilt=tilt, polygon_boundaries=polygon, - normal_rotation=axis_tilt, - reference_point=center, + roll=axis_tilt, + reference_point=reference_point if reference_point else center, ) def plot(self, ax=None, **kwargs): @@ -393,7 +516,7 @@ def plot(self, ax=None, **kwargs): Parameters ---------- - ax : matplotlib.axes.Axes, optional + 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 From 00d7ec1f99ebd1ccf9e14cf7a248b12bebb650d0 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sat, 22 Jun 2024 23:45:02 +0200 Subject: [PATCH 12/25] Hopefully last example changes --- .../plot_par_direct_shading0_fixed_tilt.py | 190 +++++++++++++----- .../plot_spatial_row_to_row_shading.py | 62 ++++-- pvlib/spatial.py | 8 +- 3 files changed, 191 insertions(+), 69 deletions(-) diff --git a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py index ae9273cc0a..366af986ce 100644 --- a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py +++ b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py @@ -52,8 +52,13 @@ import matplotlib.pyplot as plt from mpl_toolkits.mplot3d.art3d import Poly3DCollection +import matplotlib.animation as animation +from matplotlib.dates import ConciseDateFormatter, AutoDateLocator import shapely.plotting +import shapely as sp import pandas as pd +import numpy as np +from functools import partial from pvlib.spatial import RectangularSurface from pvlib import solarposition @@ -69,91 +74,176 @@ .union(fall_equinox) .union(winter_solstice) ) +# dates = spring_equinox solar_position = solarposition.get_solarposition( dates, latitude, longitude, altitude ) -solar_azimuth = solar_position["azimuth"].iloc[0] -solar_zenith = solar_position["apparent_zenith"].iloc[0] +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 top surface - center=[10, 5, 0], - azimuth=180, # chosen instead of 0 (north) for intuitive visualization +field = RectangularSurface( # crops surface + center=[5, 10, 0], + azimuth=90, tilt=0, axis_tilt=0, width=10, length=20, ) -pv_row1 = RectangularSurface( # south-most row (lowest Y-coordinate) - center=[10, -10 / 2 + 5, 2 / 2], - azimuth=180, - tilt=90, - axis_tilt=0, - width=2, - length=20, -) -pv_row2 = RectangularSurface( # north-most row (highest Y-coordinate) - center=[10, 10 / 2 + 5, 2 / 2], - azimuth=180, - tilt=90, - axis_tilt=0, - width=2, - 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, + ), ) -shades_3d = field.get_3D_shades_from( - solar_zenith, solar_azimuth, pv_row1, pv_row2 -) -shades_2d = field.get_2D_shades_from( - solar_zenith, solar_azimuth, shades_3d=shades_3d -) + +# %% +# 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 +): + print(f"{timestamp_index=}") + # 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] + # 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} - -fig = plt.figure(figsize=(10, 10)) -gs = fig.add_gridspec(10, 1) +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]) -ax1.view_init(elev=45, azim=-45) -field.plot(ax=ax1, **field_style) -pv_row1.plot(ax=ax1, **row_style) -pv_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)) - +# Upper plot, 3D ax1.axis("equal") -ax1.set_zlim(0) +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() -field_style_2d = {**field_style, "add_points": False} shapely.plotting.plot_polygon(field2d, ax=ax2, **field_style_2d) -shade_style_2d = {**shade_style, "add_points": False} -for shade in shades_2d.geoms: - shapely.plotting.plot_polygon(shade, ax=ax2, **shade_style_2d) -# %% -beam_shaded_fraction = ( - sum(shade.area for shade in shades_2d.geoms) / field2d.area +# 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( + sp.Polygon([[0, 0]] * 4), ax=ax2, **shade_style_2d + ), +) * len(pv_rows) + +ani = animation.FuncAnimation( + fig, + partial( + simulation_and_plot_callback, + shade_3d_artists=shade3d_artists, + shade_2d_artists=shade2d_artists, + ), + frames=np.arange(N), + interval=200, + blit=True, ) -print(beam_shaded_fraction) + +# %% +# Shaded Fraction vs. Time +# ------------------------ + +fig, ax = plt.subplots() + +ax.plot(dates, fixed_tilt_shaded_fraction, label="Fixed-Tilt") +locator = AutoDateLocator() +ax.xaxis.set_major_locator(locator) +ax.xaxis.set_major_formatter(ConciseDateFormatter(locator)) + +ax.set_xlabel("Time") +ax.set_ylabel("Shaded Fraction [Unitless]") + +plt.legend() +plt.show() + # %% # References diff --git a/docs/examples/shading/plot_spatial_row_to_row_shading.py b/docs/examples/shading/plot_spatial_row_to_row_shading.py index 470b58c819..cf2d5bc35c 100644 --- a/docs/examples/shading/plot_spatial_row_to_row_shading.py +++ b/docs/examples/shading/plot_spatial_row_to_row_shading.py @@ -6,12 +6,17 @@ # %% # 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 shading fraction in 3D and 2D space with the help of the +# the shaded fraction in 3D and 2D space with the help of the # :py:mod:`~pvlib.spatial` module and its classes. # -# First section of the example will show how to calculate the shading fraction -# for an instantaneous time. The second section will show how to calculate the -# shading fraction for a time range. +# 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 @@ -19,42 +24,60 @@ import shapely.plotting # %% +# 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. -solar_azimuth = 180 # degrees +solar_azimuth = 165 # degrees solar_zenith = 75 # degrees -# Create the first row of panels row1 = RectangularSurface( # south-most row center=[0, 0, 3], azimuth=165, tilt=20, axis_tilt=10, width=2, length=20 ) -# Create the second row of panels row2 = RectangularSurface( # north-most row center=[0, 3, 3], azimuth=165, tilt=20, axis_tilt=10, width=2, length=20 ) -# Calculate the shadow +# %% +# 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, row1) +shades_2d = row2.get_2D_shades_from( + solar_zenith, solar_azimuth, shades_3d=shades_3d +) # %% -# Plot the scene and the shadow -# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# 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)) -gs = fig.add_gridspec(10, 1) -ax1 = fig.add_subplot(gs[0:6, 0], projection="3d") +# 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=45, azim=-45) row1.plot(ax=ax1, **row_style) row2.plot(ax=ax1, **row_style) @@ -71,10 +94,17 @@ ax1.set_xlabel("West(-) East(+) [m]") ax1.set_ylabel("South(-) North(+) [m]") -row_style_2d = {**row_style, "add_points": False} +# 2D plot row2_2d = row2.representation_in_2D_space() -print(f"{row2_2d=}") shapely.plotting.plot_polygon(row2_2d, ax=ax2, **row_style_2d) -shade_style_2d = {**shade_style, "add_points": False} 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.area +print(f"The shaded fraction is {shaded_fraction:.2f}") diff --git a/pvlib/spatial.py b/pvlib/spatial.py index df65412bbc..d431ad8e24 100644 --- a/pvlib/spatial.py +++ b/pvlib/spatial.py @@ -351,7 +351,7 @@ def get_2D_shades_from( """ shades = self.get_3D_shades_from(solar_zenith, solar_azimuth, *others) if shades_3d is not None: - shades = sp.MultiPolygon(shades, shades_3d) + shades = sp.MultiPolygon((shades.geoms, shades_3d.geoms)) # and clip the 2D shades to the surface in 2D _self_projected_polygon = self.representation_in_2D_space() @@ -475,6 +475,7 @@ class RectangularSurface(FlatSurface): Point to use as reference for 2D projections. If not provided, the central point is used. """ + def __init__( self, center, @@ -528,5 +529,6 @@ def plot(self, ax=None, **kwargs): x, y, z = np.hsplit( np.array(self._polygon.exterior.coords[:-1]).flatten(order="F"), 3 ) - ax.plot_trisurf(x, y, z, triangles=((0, 1, 2), (0, 2, 3)), **kwargs) - return ax + return ax.plot_trisurf( + x, y, z, triangles=((0, 1, 2), (0, 2, 3)), **kwargs + ) From 7be74c83856929d5dd147414af2eadcebb1f7092 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sat, 22 Jun 2024 23:50:31 +0200 Subject: [PATCH 13/25] Update plot_spatial_row_to_row_shading.py --- docs/examples/shading/plot_spatial_row_to_row_shading.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/docs/examples/shading/plot_spatial_row_to_row_shading.py b/docs/examples/shading/plot_spatial_row_to_row_shading.py index cf2d5bc35c..e251e6b2cb 100644 --- a/docs/examples/shading/plot_spatial_row_to_row_shading.py +++ b/docs/examples/shading/plot_spatial_row_to_row_shading.py @@ -32,6 +32,8 @@ # 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 @@ -78,7 +80,10 @@ ax2 = fig.add_subplot(gs[8:, 0]) # 3D plot -ax1.view_init(elev=45, azim=-45) +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: @@ -106,5 +111,5 @@ # 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.area +shaded_fraction = sum(shade.area for shade in shades_2d.geoms) / row2_2d.area print(f"The shaded fraction is {shaded_fraction:.2f}") From 0220c29432969c2c05348aaa1913ed18e782afad Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sun, 23 Jun 2024 05:03:27 +0200 Subject: [PATCH 14/25] Improve plots --- .../plot_par_direct_shading0_fixed_tilt.py | 42 ++++++++++++++----- 1 file changed, 31 insertions(+), 11 deletions(-) diff --git a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py index 366af986ce..f79cf4d8b2 100644 --- a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py +++ b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py @@ -53,7 +53,7 @@ import matplotlib.pyplot as plt from mpl_toolkits.mplot3d.art3d import Poly3DCollection import matplotlib.animation as animation -from matplotlib.dates import ConciseDateFormatter, AutoDateLocator +from matplotlib.dates import DateFormatter import shapely.plotting import shapely as sp import pandas as pd @@ -131,12 +131,15 @@ # Shades callback def simulation_and_plot_callback( - timestamp_index, *, shade_3d_artists, shade_2d_artists + timestamp_index, *, shade_3d_artists, shade_2d_artists, sun_annotation ): - print(f"{timestamp_index=}") # 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 @@ -214,6 +217,7 @@ def simulation_and_plot_callback( sp.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, @@ -221,25 +225,41 @@ def simulation_and_plot_callback( 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, ) - +ani.save("fixed_tilt_shading.gif", writer="pillow", fps=5) # %% # Shaded Fraction vs. Time # ------------------------ +# Create a handy pandas series to plot the shaded fraction vs. time. -fig, ax = plt.subplots() +fixed_tilt_shaded_fraction = pd.Series(fixed_tilt_shaded_fraction, index=dates) -ax.plot(dates, fixed_tilt_shaded_fraction, label="Fixed-Tilt") -locator = AutoDateLocator() -ax.xaxis.set_major_locator(locator) -ax.xaxis.set_major_formatter(ConciseDateFormatter(locator)) +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) -ax.set_xlabel("Time") -ax.set_ylabel("Shaded Fraction [Unitless]") +axs[0].set_ylabel("Shaded Fraction [Unitless]") +axs[0].set_ylim(0, 1) plt.legend() plt.show() From 23e9cc22c4b0bd8a84776691b31bc5a6d4181b97 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sun, 23 Jun 2024 05:08:49 +0200 Subject: [PATCH 15/25] Update .gitignore --- .gitignore | 2 -- 1 file changed, 2 deletions(-) diff --git a/.gitignore b/.gitignore index 6465a1d3bb..1605e1780f 100644 --- a/.gitignore +++ b/.gitignore @@ -97,5 +97,3 @@ coverage.xml # airspeed velocity env results - -/.venv From 07d529d338101bcc4eff668222b047bcacc08e97 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sun, 23 Jun 2024 05:16:09 +0200 Subject: [PATCH 16/25] Update plot_par_direct_shading0_fixed_tilt.py --- .../agrivoltaics/plot_par_direct_shading0_fixed_tilt.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py index f79cf4d8b2..8fd7b670e9 100644 --- a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py +++ b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py @@ -231,7 +231,7 @@ def simulation_and_plot_callback( interval=200, blit=True, ) -ani.save("fixed_tilt_shading.gif", writer="pillow", fps=5) + # %% # Shaded Fraction vs. Time # ------------------------ From b23afb88841ddeba980e847dc22a0febcb1e17ca Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sun, 23 Jun 2024 05:43:58 +0200 Subject: [PATCH 17/25] Update test_spatial.py --- pvlib/tests/test_spatial.py | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/pvlib/tests/test_spatial.py b/pvlib/tests/test_spatial.py index 312c741717..a84a0a528a 100644 --- a/pvlib/tests/test_spatial.py +++ b/pvlib/tests/test_spatial.py @@ -131,12 +131,10 @@ def test_RectangularSurface__init__(): width=width, length=length, ) - assert surface.center == center + assert surface.reference_point == center assert surface.azimuth == surface_azimuth assert surface.tilt == surface_tilt - assert surface.axis_tilt == axis_tilt - assert surface.width == width - assert surface.length == length + assert surface.roll == axis_tilt # construct from coordinates center = [0, 0, 0] surface_azimuth = 180 @@ -152,12 +150,10 @@ def test_RectangularSurface__init__(): width=width, length=length, ) - assert surface.center == center + assert surface.reference_point == center assert surface.azimuth == surface_azimuth assert surface.tilt == surface_tilt - assert surface.axis_tilt == axis_tilt - assert surface.width == width - assert surface.length == length + assert surface.roll == axis_tilt def test_RectangularSurface__calc_surface_tilt_and_azimuth(): @@ -175,4 +171,4 @@ def test_RectangularSurface__calc_surface_tilt_and_azimuth(): length=length, ) assert surface.tilt == tilt - assert surface.azimuth == azimuth \ No newline at end of file + assert surface.azimuth == azimuth From 05398973dd688db28f9789a5dc288e479554597c Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sun, 23 Jun 2024 05:45:47 +0200 Subject: [PATCH 18/25] Add last TODO --- pvlib/spatial.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pvlib/spatial.py b/pvlib/spatial.py index d431ad8e24..f08685867b 100644 --- a/pvlib/spatial.py +++ b/pvlib/spatial.py @@ -436,7 +436,8 @@ def combine_2D_shades(self, *shades): shapely.MultiPolygon Combined shades. """ - return sp.ops.unary_union(shades) + # TODO: implement this method + return None def plot(self, ax=None, **kwargs): pass # TODO: implement this method From 7ea1a1211730f678a99ba27c476359d5cd749581 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sun, 23 Jun 2024 06:01:55 +0200 Subject: [PATCH 19/25] add version and shapely dpendency --- ci/requirements-py3.10.yml | 1 + ci/requirements-py3.11.yml | 1 + ci/requirements-py3.12.yml | 1 + ci/requirements-py3.8-min.yml | 1 + ci/requirements-py3.8.yml | 1 + ci/requirements-py3.9.yml | 1 + pyproject.toml | 2 +- 7 files changed, 7 insertions(+), 1 deletion(-) diff --git a/ci/requirements-py3.10.yml b/ci/requirements-py3.10.yml index fcb89d2e7f..2cfa224ee9 100644 --- a/ci/requirements-py3.10.yml +++ b/ci/requirements-py3.10.yml @@ -22,6 +22,7 @@ dependencies: - pytz - requests - scipy >= 1.6.0 + - shapely >= 2.0.2 - statsmodels - pip: - nrel-pysam>=2.0 diff --git a/ci/requirements-py3.11.yml b/ci/requirements-py3.11.yml index f6556ecf94..dce73e956f 100644 --- a/ci/requirements-py3.11.yml +++ b/ci/requirements-py3.11.yml @@ -22,6 +22,7 @@ dependencies: - pytz - requests - scipy >= 1.6.0 + - shapely >= 2.0.2 - statsmodels - pip: - nrel-pysam>=2.0 diff --git a/ci/requirements-py3.12.yml b/ci/requirements-py3.12.yml index f3d8fc2d0c..1f9338289f 100644 --- a/ci/requirements-py3.12.yml +++ b/ci/requirements-py3.12.yml @@ -22,6 +22,7 @@ dependencies: - pytz - requests - scipy >= 1.6.0 + - shapely >= 2.0.2 - statsmodels - pip: - nrel-pysam>=2.0 diff --git a/ci/requirements-py3.8-min.yml b/ci/requirements-py3.8-min.yml index ef230b1627..549578468c 100644 --- a/ci/requirements-py3.8-min.yml +++ b/ci/requirements-py3.8-min.yml @@ -3,6 +3,7 @@ channels: - defaults dependencies: - coveralls + - shapely >= 2.0.2 - pip - pytest - pytest-cov diff --git a/ci/requirements-py3.8.yml b/ci/requirements-py3.8.yml index 19c3b8f2f3..ce85406c92 100644 --- a/ci/requirements-py3.8.yml +++ b/ci/requirements-py3.8.yml @@ -22,6 +22,7 @@ dependencies: - pytz - requests - scipy >= 1.6.0 + - shapely >= 2.0.2 - statsmodels - pip: - nrel-pysam>=2.0 diff --git a/ci/requirements-py3.9.yml b/ci/requirements-py3.9.yml index b5aa976b4b..f03e90c584 100644 --- a/ci/requirements-py3.9.yml +++ b/ci/requirements-py3.9.yml @@ -22,6 +22,7 @@ dependencies: - pytz - requests - scipy >= 1.6.0 + - shapely >= 2.0.2 - statsmodels - pip: - nrel-pysam>=2.0 diff --git a/pyproject.toml b/pyproject.toml index 273cf2a2a7..320b36a137 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,7 +53,7 @@ optional = [ 'numba >= 0.17.0', 'solarfactors', 'statsmodels', - 'shapely', + 'shapely >= 2.0.2', ] doc = [ 'ipython', From 98d426f457c2a5a8c3e790045935d9c8cdc12322 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sun, 23 Jun 2024 23:54:53 +0200 Subject: [PATCH 20/25] Make shpaley a mandatory dependency so I can test in PR? --- ci/requirements-py3.10.yml | 1 + ci/requirements-py3.11.yml | 1 + ci/requirements-py3.12.yml | 1 + ci/requirements-py3.8-min.yml | 2 +- ci/requirements-py3.8.yml | 2 +- ci/requirements-py3.9.yml | 3 ++- 6 files changed, 7 insertions(+), 3 deletions(-) diff --git a/ci/requirements-py3.10.yml b/ci/requirements-py3.10.yml index 2cfa224ee9..89cb593ff7 100644 --- a/ci/requirements-py3.10.yml +++ b/ci/requirements-py3.10.yml @@ -27,3 +27,4 @@ dependencies: - 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 dce73e956f..5be297e6ab 100644 --- a/ci/requirements-py3.11.yml +++ b/ci/requirements-py3.11.yml @@ -27,3 +27,4 @@ dependencies: - 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 1f9338289f..489350aac0 100644 --- a/ci/requirements-py3.12.yml +++ b/ci/requirements-py3.12.yml @@ -27,3 +27,4 @@ dependencies: - 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 549578468c..aed65ddb32 100644 --- a/ci/requirements-py3.8-min.yml +++ b/ci/requirements-py3.8-min.yml @@ -3,7 +3,6 @@ channels: - defaults dependencies: - coveralls - - shapely >= 2.0.2 - pip - pytest - pytest-cov @@ -20,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 ce85406c92..0ed6bfa119 100644 --- a/ci/requirements-py3.8.yml +++ b/ci/requirements-py3.8.yml @@ -22,8 +22,8 @@ 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.9.yml b/ci/requirements-py3.9.yml index f03e90c584..19a7194b7d 100644 --- a/ci/requirements-py3.9.yml +++ b/ci/requirements-py3.9.yml @@ -26,4 +26,5 @@ dependencies: - statsmodels - pip: - nrel-pysam>=2.0 - - solarfactors \ No newline at end of file + - solarfactors + - shapely >= 2.0.2 From dbc1aac056883fbd41938f79f3aa57d853f1768a Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Sun, 23 Jun 2024 23:55:43 +0200 Subject: [PATCH 21/25] Minor update to plot --- .../agrivoltaics/plot_par_direct_shading0_fixed_tilt.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py index 8fd7b670e9..f25d02c9ec 100644 --- a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py +++ b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py @@ -232,6 +232,11 @@ def simulation_and_plot_callback( blit=True, ) +# uncomment to run and save animation locally +# ani.save("fixed_tilt_shading.gif", writer="pillow") + +plt.show() + # %% # Shaded Fraction vs. Time # ------------------------ @@ -261,7 +266,6 @@ def simulation_and_plot_callback( axs[0].set_ylabel("Shaded Fraction [Unitless]") axs[0].set_ylim(0, 1) -plt.legend() plt.show() From 927aec8104d58e44d893f8de09fbf9c543ca9fb3 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Thu, 27 Jun 2024 23:43:55 +0200 Subject: [PATCH 22/25] FINALLY SHAPELY 2/1.8.5 INDEPENDENT CODE HAHAHAHAHA --- .../plot_par_direct_shading0_fixed_tilt.py | 9 +- .../plot_spatial_row_to_row_shading.py | 2 +- pvlib/spatial.py | 219 ++++++++++++++++-- pvlib/tests/test_spatial.py | 12 +- pyproject.toml | 1 + 5 files changed, 216 insertions(+), 27 deletions(-) diff --git a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py index f25d02c9ec..c7b1593ccd 100644 --- a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py +++ b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py @@ -54,8 +54,11 @@ from mpl_toolkits.mplot3d.art3d import Poly3DCollection import matplotlib.animation as animation from matplotlib.dates import DateFormatter -import shapely.plotting -import shapely as sp +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 @@ -214,7 +217,7 @@ def simulation_and_plot_callback( ) * len(pv_rows) shade2d_artists = ( shapely.plotting.plot_polygon( - sp.Polygon([[0, 0]] * 4), ax=ax2, **shade_style_2d + 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") diff --git a/docs/examples/shading/plot_spatial_row_to_row_shading.py b/docs/examples/shading/plot_spatial_row_to_row_shading.py index e251e6b2cb..a765688ca6 100644 --- a/docs/examples/shading/plot_spatial_row_to_row_shading.py +++ b/docs/examples/shading/plot_spatial_row_to_row_shading.py @@ -21,7 +21,7 @@ from pvlib.spatial import RectangularSurface import matplotlib.pyplot as plt from mpl_toolkits.mplot3d.art3d import Poly3DCollection -import shapely.plotting +import shapely # %% # Description of the system diff --git a/pvlib/spatial.py b/pvlib/spatial.py index f08685867b..08f6e9e95c 100644 --- a/pvlib/spatial.py +++ b/pvlib/spatial.py @@ -4,10 +4,196 @@ from pvlib.tools import sind, cosd import numpy as np -import shapely as sp 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): """ @@ -169,7 +355,7 @@ def __init__( self._tilt = np.array(tilt) self._roll = np.array(roll) # works for polygon_boundaries := array[N, 3] | shapely.Polygon - self._polygon = sp.Polygon(polygon_boundaries) + self._polygon = Polygon(polygon_boundaries) # internal 2D coordinates-system to translate projections matrix # only defined if needed later on self._reference_point = ( @@ -291,12 +477,15 @@ def get_3D_shade_from_flat_surface(other): count=len(coords_to_project), # speeds up allocation ) # create shapely shade object and bound it to the surface - shade = sp.Polygon(projected_vertices).intersection( - _polygon, grid_size=1e-12 - ) - return shade if isinstance(shade, sp.Polygon) else sp.Polygon() + 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 sp.MultiPolygon(map(get_3D_shade_from_flat_surface, others)) + return MultiPolygon(map(get_3D_shade_from_flat_surface, others)) def get_2D_shades_from( self, solar_zenith, solar_azimuth, *others, shades_3d=None @@ -351,7 +540,7 @@ def get_2D_shades_from( """ shades = self.get_3D_shades_from(solar_zenith, solar_azimuth, *others) if shades_3d is not None: - shades = sp.MultiPolygon((shades.geoms, shades_3d.geoms)) + 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() @@ -361,15 +550,13 @@ def get_2D_shade_from_flat_surface(other): # create shapely shade object and bound it to the surface projected_vertices_2d = self._transform_to_2D(coords_to_project) - shade = sp.Polygon(projected_vertices_2d).intersection( + shade = Polygon(projected_vertices_2d).intersection( _self_projected_polygon ) # Return Polygons only, geometries with empty area are omitted - return shade if isinstance(shade, sp.Polygon) else sp.Polygon() + return shade if isinstance(shade, Polygon) else Polygon() - return sp.MultiPolygon( - map(get_2D_shade_from_flat_surface, shades.geoms) - ) + return MultiPolygon(map(get_2D_shade_from_flat_surface, shades.geoms)) def _get_projection_transform(self): """ @@ -416,9 +603,7 @@ def representation_in_2D_space(self): """ if self._projected_polygon is None: vertices = np.array(self._polygon.exterior.coords[:-1]) - self._projected_polygon = sp.Polygon( - self._transform_to_2D(vertices) - ) + self._projected_polygon = Polygon(self._transform_to_2D(vertices)) return self._projected_polygon def combine_2D_shades(self, *shades): @@ -503,7 +688,7 @@ def __init__( _rotation = Rotation.from_euler( "ZXZ", [-azimuth, -tilt, axis_tilt], degrees=True ) - polygon = sp.Polygon(_rotation.apply(corners) + center) + polygon = Polygon(_rotation.apply(corners) + center) super().__init__( azimuth=azimuth, tilt=tilt, diff --git a/pvlib/tests/test_spatial.py b/pvlib/tests/test_spatial.py index a84a0a528a..7ddcd5eed0 100644 --- a/pvlib/tests/test_spatial.py +++ b/pvlib/tests/test_spatial.py @@ -2,7 +2,7 @@ from pvlib import spatial from numpy.testing import assert_allclose -import shapely +from shapely.geometry import Polygon import pytest @@ -59,7 +59,7 @@ def test_FlatSurface__init__(): # construct with native shapely polygon surface_azimuth = 180 surface_tilt = 0 - polygon = shapely.Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) + polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) surface = spatial.FlatSurface( azimuth=surface_azimuth, tilt=surface_tilt, @@ -68,7 +68,7 @@ def test_FlatSurface__init__(): assert surface.azimuth == surface_azimuth assert surface.tilt == surface_tilt assert surface.polygon == polygon - assert isinstance(surface.polygon, shapely.Polygon) + assert isinstance(surface.polygon, Polygon) # construct from coordinates surface_azimuth = 180 surface_tilt = 0 @@ -78,14 +78,14 @@ def test_FlatSurface__init__(): ) assert surface.azimuth == surface_azimuth assert surface.tilt == surface_tilt - assert surface.polygon == shapely.Polygon(polygon) - assert isinstance(surface.polygon, shapely.Polygon) + assert surface.polygon == Polygon(polygon) + assert isinstance(surface.polygon, Polygon) def test_FlatSurface_readonly_properties(): surface_azimuth = 180 surface_tilt = 0 - polygon = shapely.Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) + polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]) surface = spatial.FlatSurface( azimuth=surface_azimuth, tilt=surface_tilt, diff --git a/pyproject.toml b/pyproject.toml index 320b36a137..a96d2883ad 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -65,6 +65,7 @@ doc = [ 'pillow', 'sphinx-toggleprompt >= 0.0.5', 'solarfactors', + 'shapely', ] test = [ 'pytest', From ed8059038559971f6d43f5b3063bea776fba3907 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Fri, 28 Jun 2024 00:38:09 +0200 Subject: [PATCH 23/25] Be consistent with problem description --- docs/examples/shading/plot_spatial_row_to_row_shading.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/examples/shading/plot_spatial_row_to_row_shading.py b/docs/examples/shading/plot_spatial_row_to_row_shading.py index a765688ca6..9fa77d8079 100644 --- a/docs/examples/shading/plot_spatial_row_to_row_shading.py +++ b/docs/examples/shading/plot_spatial_row_to_row_shading.py @@ -43,7 +43,7 @@ ) row2 = RectangularSurface( # north-most row - center=[0, 3, 3], azimuth=165, tilt=20, axis_tilt=10, width=2, length=20 + center=[0, 3, 3], azimuth=165, tilt=30, axis_tilt=10, width=2, length=20 ) # %% From d6f38661a9f0cd0e94116bc4c0b3641f0faac978 Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Fri, 28 Jun 2024 00:41:09 +0200 Subject: [PATCH 24/25] Escape bullet points in table --- .../agrivoltaics/plot_par_direct_shading0_fixed_tilt.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py index c7b1593ccd..c69b1fc8b4 100644 --- a/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py +++ b/docs/examples/agrivoltaics/plot_par_direct_shading0_fixed_tilt.py @@ -35,17 +35,17 @@ # +--------------------+----------+-------------+----------+-------+ # | Crop area | 200 | 200 | 200 | [m²] | # +--------------------+----------+-------------+----------+-------+ -# | Pitch | - | - | 2 | [m] | +# | Pitch | \- | \- | 2 | [m] | # +--------------------+----------+-------------+----------+-------+ # | Height | 0 | 3 | 3 | [m] | # +--------------------+----------+-------------+----------+-------+ -# | Fixed tilt angle | 90 | - | - | [°] | +# | Fixed tilt angle | 90 | \- | \- | [°] | # +--------------------+----------+-------------+----------+-------+ # | Azimuth angle | 0 | 0 | 0 | [°] | # +--------------------+----------+-------------+----------+-------+ -# | Maximum tilt angle | - | 60 | 60 | [°] | +# | Maximum tilt angle | \- | 60 | 60 | [°] | # +--------------------+----------+-------------+----------+-------+ -# | Minimum tilt angle | - | -60 | -60 | [°] | +# | Minimum tilt angle | \- | -60 | -60 | [°] | # +--------------------+----------+-------------+----------+-------+ # # From 520b33e52f403ce59b6d461351638dd7e8e9749c Mon Sep 17 00:00:00 2001 From: echedey-ls <80125792+echedey-ls@users.noreply.github.com> Date: Fri, 28 Jun 2024 00:47:33 +0200 Subject: [PATCH 25/25] Add items to docs --- docs/sphinx/source/reference/index.rst | 1 + docs/sphinx/source/reference/spatial.rst | 12 ++++++++++++ 2 files changed, 13 insertions(+) create mode 100644 docs/sphinx/source/reference/spatial.rst 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