"""
The ``shading`` module contains functions that model module shading and the
associated effects on PV module output
"""
import numpy as np
import pandas as pd
from pvlib.tools import sind, cosd
[docs]
def ground_angle(surface_tilt, gcr, slant_height):
"""
Angle from horizontal of the line from a point on the row slant length
to the bottom of the facing row.
The angles are clockwise from horizontal, rather than the usual
counterclockwise direction.
Parameters
----------
surface_tilt : numeric
Surface tilt angle in degrees from horizontal, e.g., surface facing up
= 0, surface facing horizon = 90. [degree]
gcr : float
ground coverage ratio, ratio of row slant length to row spacing.
[unitless]
slant_height : numeric
The distance up the module's slant height to evaluate the ground
angle, as a fraction [0-1] of the module slant height [unitless].
Returns
-------
psi : numeric
Angle [degree].
"""
# : \\ \
# : \\ \
# : \\ \
# : \\ \ facing row
# : \\.___________\
# : \ ^*-. psi \
# : \ x *-. \
# : \ v *-.\
# : \<-----P---->\
x1 = gcr * slant_height * sind(surface_tilt)
x2 = gcr * slant_height * cosd(surface_tilt) + 1
psi = np.arctan2(x1, x2) # do this before rad2deg because it handles 0 / 0
return np.rad2deg(psi)
[docs]
def masking_angle(surface_tilt, gcr, slant_height):
"""
The elevation angle below which diffuse irradiance is blocked.
The ``height`` parameter determines how far up the module's surface to
evaluate the masking angle. The lower the point, the steeper the masking
angle [1]_. SAM uses a "worst-case" approach where the masking angle
is calculated for the bottom of the array (i.e. ``slant_height=0``) [2]_.
Parameters
----------
surface_tilt : numeric
Panel tilt from horizontal [degrees].
gcr : float
The ground coverage ratio of the array [unitless].
slant_height : numeric
The distance up the module's slant height to evaluate the masking
angle, as a fraction [0-1] of the module slant height [unitless].
Returns
-------
mask_angle : numeric
Angle from horizontal where diffuse light is blocked by the
preceding row [degrees].
See Also
--------
masking_angle_passias
sky_diffuse_passias
References
----------
.. [1] D. Passias and B. Källbäck, "Shading effects in rows of solar cell
panels", Solar Cells, Volume 11, Pages 281-291. 1984.
:doi:`10.1016/0379-6787(84)90017-6`
.. [2] Gilman, P. et al., (2018). "SAM Photovoltaic Model Technical
Reference Update", NREL Technical Report NREL/TP-6A20-67399.
Available at https://www.nrel.gov/docs/fy18osti/67399.pdf
"""
# The original equation (8 in [1]) requires pitch and collector width,
# but it's easy to non-dimensionalize it to make it a function of GCR
# by factoring out B from the argument to arctan.
numerator = gcr * (1 - slant_height) * sind(surface_tilt)
denominator = 1 - gcr * (1 - slant_height) * cosd(surface_tilt)
phi = np.arctan(numerator / denominator)
return np.degrees(phi)
[docs]
def masking_angle_passias(surface_tilt, gcr):
r"""
The average masking angle over the slant height of a row.
The masking angle is the angle from horizontal where the sky dome is
blocked by the row in front. The masking angle is larger near the lower
edge of a row than near the upper edge. This function calculates the
average masking angle as described in [1]_.
Parameters
----------
surface_tilt : numeric
Panel tilt from horizontal [degrees].
gcr : float
The ground coverage ratio of the array [unitless].
Returns
----------
mask_angle : numeric
Average angle from horizontal where diffuse light is blocked by the
preceding row [degrees].
See Also
--------
masking_angle
sky_diffuse_passias
Notes
-----
The pvlib-python authors believe that Eqn. 9 in [1]_ is incorrect.
Here we use an independent equation. First, Eqn. 8 is non-dimensionalized
(recasting in terms of GCR):
.. math::
\psi(z') = \arctan \left [
\frac{(1 - z') \sin \beta}
{\mathrm{GCR}^{-1} + (z' - 1) \cos \beta}
\right ]
Where :math:`GCR = B/C` and :math:`z' = z/B`. The average masking angle
:math:`\overline{\psi} = \int_0^1 \psi(z') \mathrm{d}z'` is then
evaluated symbolically using Maxima (using :math:`X = 1/\mathrm{GCR}`):
.. code-block:: none
load(scifac) /* for the gcfac function */
assume(X>0, cos(beta)>0, cos(beta)-X<0); /* X is 1/GCR */
gcfac(integrate(atan((1-z)*sin(beta)/(X+(z-1)*cos(beta))), z, 0, 1))
This yields the equation implemented by this function:
.. math::
\overline{\psi} = \
&-\frac{X}{2} \sin\beta \log | 2 X \cos\beta - (X^2 + 1)| \\
&+ (X \cos\beta - 1) \arctan \frac{X \cos\beta - 1}{X \sin\beta} \\
&+ (1 - X \cos\beta) \arctan \frac{\cos\beta}{\sin\beta} \\
&+ X \log X \sin\beta
The pvlib-python authors have validated this equation against numerical
integration of :math:`\overline{\psi} = \int_0^1 \psi(z') \mathrm{d}z'`.
References
----------
.. [1] D. Passias and B. Källbäck, "Shading effects in rows of solar cell
panels", Solar Cells, Volume 11, Pages 281-291. 1984.
:doi:`10.1016/0379-6787(84)90017-6`
"""
# wrap it in an array so that division by zero is handled well
beta = np.radians(np.array(surface_tilt))
sin_b = np.sin(beta)
cos_b = np.cos(beta)
X = 1/gcr
with np.errstate(divide='ignore', invalid='ignore'): # ignore beta=0
term1 = -X * sin_b * np.log(np.abs(2 * X * cos_b - (X**2 + 1))) / 2
term2 = (X * cos_b - 1) * np.arctan((X * cos_b - 1) / (X * sin_b))
term3 = (1 - X * cos_b) * np.arctan(cos_b / sin_b)
term4 = X * np.log(X) * sin_b
psi_avg = term1 + term2 + term3 + term4
# when beta=0, divide by zero makes psi_avg NaN. replace with 0:
psi_avg = np.where(np.isfinite(psi_avg), psi_avg, 0)
if isinstance(surface_tilt, pd.Series):
psi_avg = pd.Series(psi_avg, index=surface_tilt.index)
return np.degrees(psi_avg)
[docs]
def sky_diffuse_passias(masking_angle):
r"""
The diffuse irradiance loss caused by row-to-row sky diffuse shading.
Even when the sun is high in the sky, a row's view of the sky dome will
be partially blocked by the row in front. This causes a reduction in the
diffuse irradiance incident on the module. The reduction depends on the
masking angle, the elevation angle from a point on the shaded module to
the top of the shading row. In [1]_ the masking angle is calculated as
the average across the module height. SAM assumes the "worst-case" loss
where the masking angle is calculated for the bottom of the array [2]_.
This function, as in [1]_, makes the assumption that sky diffuse
irradiance is isotropic.
Parameters
----------
masking_angle : numeric
The elevation angle below which diffuse irradiance is blocked
[degrees].
Returns
-------
derate : numeric
The fraction [0-1] of blocked sky diffuse irradiance.
See Also
--------
masking_angle
masking_angle_passias
References
----------
.. [1] D. Passias and B. Källbäck, "Shading effects in rows of solar cell
panels", Solar Cells, Volume 11, Pages 281-291. 1984.
:doi:`10.1016/0379-6787(84)90017-6`
.. [2] Gilman, P. et al., (2018). "SAM Photovoltaic Model Technical
Reference Update", NREL Technical Report NREL/TP-6A20-67399.
Available at https://www.nrel.gov/docs/fy18osti/67399.pdf
"""
return 1 - cosd(masking_angle/2)**2
[docs]
def projected_solar_zenith_angle(solar_zenith, solar_azimuth,
axis_tilt, axis_azimuth):
r"""
Calculate projected solar zenith angle in degrees.
This solar zenith angle is projected onto the plane whose normal vector is
defined by ``axis_tilt`` and ``axis_azimuth``. The normal vector is in the
direction of ``axis_azimuth`` (clockwise from north) and tilted from
horizontal by ``axis_tilt``. See Figure 5 in [1]_:
.. figure:: ../../_images/Anderson_Mikofski_2020_Fig5.jpg
:alt: Wire diagram of coordinates systems to obtain the projected angle.
:align: center
:scale: 50 %
Fig. 5, [1]_: Solar coordinates projection onto tracker rotation plane.
Parameters
----------
solar_zenith : numeric
Sun's apparent zenith in degrees.
solar_azimuth : numeric
Sun's azimuth in degrees.
axis_tilt : numeric
Axis tilt angle in degrees. From horizontal plane to array plane.
axis_azimuth : numeric
Axis azimuth angle in degrees.
North = 0°; East = 90°; South = 180°; West = 270°
Returns
-------
Projected_solar_zenith : numeric
In degrees.
Notes
-----
This projection has a variety of applications in PV. For example:
- Projecting the sun's position onto the plane perpendicular to
the axis of a single-axis tracker (i.e. the plane
whose normal vector coincides with the tracker torque tube)
yields the tracker rotation angle that maximizes direct irradiance
capture. This tracking strategy is called *true-tracking*. Learn more
about tracking in
:ref:`sphx_glr_gallery_solar-tracking_plot_single_axis_tracking.py`.
- Self-shading in large PV arrays is often modeled by assuming
a simplified 2-D array geometry where the sun's position is
projected onto the plane perpendicular to the PV rows.
The projected zenith angle is then used for calculations
regarding row-to-row shading.
Examples
--------
Calculate the ideal true-tracking angle for a horizontal north-south
single-axis tracker:
>>> rotation = projected_solar_zenith_angle(solar_zenith, solar_azimuth,
>>> axis_tilt=0, axis_azimuth=180)
Calculate the projected zenith angle in a south-facing fixed tilt array
(note: the ``axis_azimuth`` of a fixed-tilt row points along the length
of the row):
>>> psza = projected_solar_zenith_angle(solar_zenith, solar_azimuth,
>>> axis_tilt=0, axis_azimuth=90)
References
----------
.. [1] K. Anderson and M. Mikofski, 'Slope-Aware Backtracking for
Single-Axis Trackers', National Renewable Energy Lab. (NREL), Golden,
CO (United States);
NREL/TP-5K00-76626, Jul. 2020. :doi:`10.2172/1660126`.
See Also
--------
pvlib.solarposition.get_solarposition
"""
# Assume the tracker reference frame is right-handed. Positive y-axis is
# oriented along tracking axis; from north, the y-axis is rotated clockwise
# by the axis azimuth and tilted from horizontal by the axis tilt. The
# positive x-axis is 90 deg clockwise from the y-axis and parallel to
# horizontal (e.g., if the y-axis is south, the x-axis is west); the
# positive z-axis is normal to the x and y axes, pointed upward.
# Since elevation = 90 - zenith, sin(90-x) = cos(x) & cos(90-x) = sin(x):
# Notation from [1], modified to use zenith instead of elevation
# cos(elevation) = sin(zenith) and sin(elevation) = cos(zenith)
# Avoid recalculating these values
sind_solar_zenith = sind(solar_zenith)
cosd_axis_azimuth = cosd(axis_azimuth)
sind_axis_azimuth = sind(axis_azimuth)
sind_axis_tilt = sind(axis_tilt)
# Sun's x, y, z coords
sx = sind_solar_zenith * sind(solar_azimuth)
sy = sind_solar_zenith * cosd(solar_azimuth)
sz = cosd(solar_zenith)
# Eq. (4); sx', sz' values from sun coordinates projected onto surface
sx_prime = sx * cosd_axis_azimuth - sy * sind_axis_azimuth
sz_prime = (
sx * sind_axis_azimuth * sind_axis_tilt
+ sy * sind_axis_tilt * cosd_axis_azimuth
+ sz * cosd(axis_tilt)
)
# Eq. (5); angle between sun's beam and surface
theta_T = np.degrees(np.arctan2(sx_prime, sz_prime))
return theta_T
[docs]
def shaded_fraction1d(
solar_zenith,
solar_azimuth,
axis_azimuth,
shaded_row_rotation,
*,
collector_width,
pitch,
axis_tilt=0,
surface_to_axis_offset=0,
cross_axis_slope=0,
shading_row_rotation=None,
):
r"""
Shaded fraction in the vertical dimension of tilted rows, or perpendicular
to the axis of horizontal rows.
If ``shading_row_rotation`` isn't provided, it is assumed that
both the shaded row and the shading row (the one blocking the
direct beam) have the same rotation and azimuth values.
.. warning::
The function assumes that the roles of the shaded and shading rows
remain the same during the day. In the case where the shading and
shaded rows change throughout the day, e.g. a N-S single-axis tracker,
the inputs must be switched depending on the sign of the projected
solar zenith angle. See the Examples section below.
.. versionadded:: 0.11.0
Parameters
----------
solar_zenith : numeric
Solar zenith angle, in degrees.
solar_azimuth : numeric
Solar azimuth angle, in degrees.
axis_azimuth : numeric
Axis azimuth of the rotation axis of a tracker, in degrees.
Fixed-tilt arrays can be considered as a particular case of a tracker.
North=0º, South=180º, East=90º, West=270º.
shaded_row_rotation : numeric
Right-handed rotation of the row receiving the shade, with respect
to ``axis_azimuth``. In degrees :math:`^{\circ}`.
collector_width : numeric
Vertical length of a tilted row. The returned ``shaded_fraction``
is the ratio of the shadow over this value.
pitch : numeric
Axis-to-axis horizontal spacing of the row.
axis_tilt : numeric, default 0
Tilt of the rows axis from horizontal. In degrees :math:`^{\circ}`.
surface_to_axis_offset : numeric, default 0
Distance between the rotating axis and the collector surface.
May be used to account for a torque tube offset.
cross_axis_slope : numeric, default 0
Angle of the plane containing the rows' axes from
horizontal. Right-handed rotation with respect to the rows axes.
In degrees :math:`^{\circ}`.
shading_row_rotation : numeric, optional
Right-handed rotation of the row casting the shadow, with respect
to the row axis. In degrees :math:`^{\circ}`.
Returns
-------
shaded_fraction : numeric
The fraction of the collector width shaded by an adjacent row. A
value of 1 is completely shaded and 0 is no shade.
Notes
-----
All length parameters must have the same units.
Parameters are defined as follow:
.. figure:: ../../_images/Anderson_Jensen_2024_Fig3.png
:alt: Diagram showing the two rows and the parameters of the model.
Figure 3 of [1]_. See correspondence between this nomenclature and the
function parameters in the table below.
+------------------+----------------------------+---------------------+
| Symbol | Parameter | Units |
+==================+============================+=====================+
| :math:`\theta_1` | ``shading_row_rotation`` | |
+------------------+----------------------------+ |
| :math:`\theta_2` | ``shaded_row_rotation`` | Degrees |
+------------------+----------------------------+ :math:`^{\circ}` |
| :math:`\beta_c` | ``cross_axis_slope`` | |
+------------------+----------------------------+---------------------+
| :math:`p` | ``pitch`` | Any consistent |
+------------------+----------------------------+ length unit across |
| :math:`\ell` | ``collector_width`` | all these |
+------------------+----------------------------+ parameters, e.g. |
| :math:`z_0` | ``surface_to_axis_offset`` | :math:`m`. |
+------------------+----------------------------+---------------------+
| :math:`f_s` | Return value | Dimensionless |
+------------------+----------------------------+---------------------+
Examples
--------
**Fixed-tilt south-facing array on flat terrain**
Tilted row with a pitch of 3 m, a collector width of
2 m, and row rotations of 30°. In the morning.
>>> shaded_fraction1d(solar_zenith=80, solar_azimuth=135,
... axis_azimuth=90, shaded_row_rotation=30, shading_row_rotation=30,
... collector_width=2, pitch=3, axis_tilt=0,
... surface_to_axis_offset=0.05, cross_axis_slope=0)
0.47755694708090535
**Fixed-tilt north-facing array on sloped terrain**
Tilted row with a pitch of 4 m, a collector width of
2.5 m, and row rotations of 50° for the shaded
row and 30° for the shading row. The rows are on a
10° slope, where their axis is on the most inclined
direction (zero cross-axis slope). Shaded in the morning.
>>> shaded_fraction1d(solar_zenith=80, solar_azimuth=75.5,
... axis_azimuth=270, shaded_row_rotation=50, shading_row_rotation=30,
... collector_width=2.5, pitch=4, axis_tilt=10,
... surface_to_axis_offset=0.05, cross_axis_slope=0)
0.793244836197256
**N-S single-axis tracker on sloped terrain**
Horizontal trackers with a pitch of 3 m, a collector width of
1.4 m, and tracker rotations of 30° pointing east,
in the morning. Terrain slope is 7° west-east (east-most
tracker is higher than the west-most tracker).
>>> shaded_fraction1d(solar_zenith=80, solar_azimuth=90, axis_azimuth=180,
... shaded_row_rotation=-30, collector_width=1.4, pitch=3, axis_tilt=0,
... surface_to_axis_offset=0.10, cross_axis_slope=7)
0.8242176864434579
Note the previous example only is valid for the shaded fraction of the
west-most tracker in the morning, and assuming it is the
shaded tracker during all the day is incorrect.
During the afternoon, it is the one casting the shadow onto the
east-most tracker.
To calculate the shaded fraction for the east-most
tracker, you must input the corresponding ``shaded_row_rotation``
in the afternoon.
>>> shaded_fraction1d(solar_zenith=80, solar_azimuth=270, axis_azimuth=180,
... shaded_row_rotation=30, collector_width=1.4, pitch=3, axis_tilt=0,
... surface_to_axis_offset=0.10, cross_axis_slope=7)
0.018002567182254348
You must switch the input/output depending on the
sign of the projected solar zenith angle. See
:py:func:`~pvlib.shading.projected_solar_zenith_angle` and the example
:ref:`sphx_glr_gallery_shading_plot_shaded_fraction1d_ns_hsat_example.py`
See also
--------
pvlib.shading.projected_solar_zenith_angle
References
----------
.. [1] 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`
"""
# For nomenclature you may refer to [1].
# rotation of row casting the shadow defaults to shaded row's one
if shading_row_rotation is None:
shading_row_rotation = shaded_row_rotation
# projected solar zenith angle
projected_solar_zenith = projected_solar_zenith_angle(
solar_zenith,
solar_azimuth,
axis_tilt,
axis_azimuth,
)
# calculate repeated elements
thetas_1_S_diff = shading_row_rotation - projected_solar_zenith
thetas_2_S_diff = shaded_row_rotation - projected_solar_zenith
thetaS_rotation_diff = projected_solar_zenith - cross_axis_slope
cos_theta_2_S_diff_abs = np.abs(cosd(thetas_2_S_diff))
# Eq. (12) of [1]
t_asterisk = (
0.5
+ np.abs(cosd(thetas_1_S_diff)) / cos_theta_2_S_diff_abs / 2
+ (
np.sign(projected_solar_zenith)
* surface_to_axis_offset
/ collector_width
/ cos_theta_2_S_diff_abs
* (sind(thetas_2_S_diff) - sind(thetas_1_S_diff))
)
- (
pitch
/ collector_width
* cosd(thetaS_rotation_diff)
/ cos_theta_2_S_diff_abs
/ cosd(cross_axis_slope)
)
)
return np.clip(t_asterisk, 0, 1)
[docs]
def direct_martinez(
poa_global,
poa_direct,
shaded_fraction,
shaded_blocks,
total_blocks,
):
r"""
A shading loss power factor for non-monolithic silicon
modules and arrays with an arbitrary number of bypass diodes.
This experimental model reduces the direct and circumsolar
irradiance reaching the module's cells based on the number of *blocks*
affected by the shadow.
More on blocks in the *Notes* section and in [1]_.
.. versionadded:: 0.11.0
Parameters
----------
poa_global : numeric
Plane of array global irradiance. [W/m²].
poa_direct : numeric
Plane of array direct and circumsolar irradiance. [W/m²].
shaded_fraction : numeric
Fraction of module surface area that is shaded. [Unitless].
shaded_blocks : numeric
Number of blocks affected by the shadow. [Unitless].
If a floating point number is provided, it will be rounded up.
total_blocks : int
Number of total blocks. Unitless.
Returns
-------
shading_losses : numeric
Fraction of DC power lost due to shading. [Unitless]
Notes
-----
The implemented equations are (6) and (8) from [1]_:
.. math::
(1 - F_{ES}) = (1 - F_{GS}) \left(1 - \frac{N_{SB}}{N_{TB} + 1}\right)
\quad \text{(6)}
\left(1 - \frac{P_{S}}{P_{NS}}\right) = \left(1 -
\frac{\left[(B + D^{CIR})(1 - F_{ES}) + D^{ISO} + R\right]}{G}\right)
\quad \text{(8)}
In (6), :math:`(1 - F_{ES})` is the correction factor to be multiplied by
the direct and circumsolar irradiance, :math:`F_{GS}` is the shaded
fraction of the collector, :math:`N_{SB}` is the number of shaded blocks
and :math:`N_{TB}` is the number of total blocks.
In (8), :math:`\frac{P_{S}}{P_{NS}}` is the fraction of DC power lost due
to shading, :math:`P_{S}` is the power output of the shaded module,
:math:`P_{NS}` is the power output of the non-shaded module,
:math:`B + D^{CIR}` is the beam and circumsolar irradiance,
:math:`D^{ISO} + R` is the sum of diffuse and albedo irradiances and
:math:`G` is the global irradiance.
**Blocks terminology:**
A *block* is defined in [1]_ as a group of solar cells protected by a
bypass diode. Also, a *block* is shaded when at least one of its
cells is partially shaded.
The total number of blocks and their layout depend on the module(s) used.
Many manufacturers don't specify this information explicitly.
However, these values can be inferred from:
- the number of bypass diodes
- where and how many junction boxes are present on the back of the module
- whether or not the module is comprised of *half-cut cells*
The latter two are heavily correlated.
For example:
1. A module with 1 bypass diode behaves as 1 block.
2. A module with 3 bypass diodes and 1 junction box is likely to have 3
blocks.
3. A half-cut module with 3 junction boxes (split junction boxes) is
likely to have 3x2 blocks. The number of blocks along the longest
side of the module is 2 and along the shortest side is 3.
4. A module without bypass diodes doesn't constitute a block, but may be
part of one.
Examples
--------
Minimal example. For a complete example, see
:ref:`sphx_glr_gallery_shading_plot_martinez_shade_loss.py`.
>>> import numpy as np
>>> from pvlib import shading
>>> total_blocks = 3 # blocks along the vertical of the module
>>> POA_direct_and_circumsolar, POA_diffuse = 600, 80 # W/m²
>>> POA_global = POA_direct_and_circumsolar + POA_diffuse
>>> P_out_unshaded = 3000 # W
>>> # calculation of the shaded fraction for the collector
>>> shaded_fraction = shading.shaded_fraction1d(
>>> solar_zenith=80, solar_azimuth=180,
>>> axis_azimuth=90, shaded_row_rotation=25,
>>> collector_width=0.5, pitch=1, surface_to_axis_offset=0,
>>> cross_axis_slope=5.711, shading_row_rotation=50)
>>> # calculation of the number of shaded blocks
>>> shaded_blocks = np.ceil(total_blocks*shaded_fraction)
>>> # apply the Martinez power losses to the calculated shading
>>> loss_fraction = shading.direct_martinez(
>>> POA_global, POA_direct_and_circumsolar,
>>> shaded_fraction, shaded_blocks, total_blocks)
>>> P_out_corrected = P_out_unshaded * (1 - loss_fraction)
See Also
--------
shaded_fraction1d : to calculate 1-dimensional shaded fraction
References
----------
.. [1] F. Martínez-Moreno, J. Muñoz, and E. Lorenzo, 'Experimental model
to estimate shading losses on PV arrays', Solar Energy Materials and
Solar Cells, vol. 94, no. 12, pp. 2298-2303, Dec. 2010,
:doi:`10.1016/j.solmat.2010.07.029`.
""" # Contributed by Echedey Luis, 2024
beam_factor = ( # Eq. (6) of [1]
(1 - shaded_fraction)
* (1 - np.ceil(shaded_blocks) / (1 + total_blocks))
)
return ( # Eq. (8) of [1]
1
- (
poa_direct * beam_factor
+ (poa_global - poa_direct) # diffuse and albedo
)
/ poa_global
)
