diff --git a/pvlib/bifacial/__init__.py b/pvlib/bifacial/__init__.py index e166c55108..256a93b99f 100644 --- a/pvlib/bifacial/__init__.py +++ b/pvlib/bifacial/__init__.py @@ -3,7 +3,9 @@ """ from pvlib._deprecation import deprecated -from pvlib.bifacial import pvfactors, infinite_sheds, utils # noqa: F401 +from pvlib.bifacial import ( # noqa: F401 + ants2d, infinite_sheds, pvfactors, utils +) from .loss_models import power_mismatch_deline # noqa: F401 pvfactors_timeseries = deprecated( diff --git a/pvlib/bifacial/ants2d.py b/pvlib/bifacial/ants2d.py new file mode 100644 index 0000000000..d06e181e98 --- /dev/null +++ b/pvlib/bifacial/ants2d.py @@ -0,0 +1,591 @@ +r""" +Functions for the ANTS 2D bifacial irradiance model. +""" + +import numpy as np +import pandas as pd +from pvlib.tools import cosd, sind, tand, acosd +from pvlib.bifacial import utils +from pvlib.irradiance import aoi_projection, haydavies, perez +from pvlib.shading import projected_solar_zenith_angle, shaded_fraction1d +from pvlib.tracking import calc_surface_orientation + + +def _shaded_fraction(tracker_rotation, phi, gcr, x0=0, x1=1): + """ + Calculate fraction (from the bottom) of row slant height that is shaded + from direct irradiance by the row in front toward the sun. + + Parameters + ---------- + tracker_rotation : numeric + Tracker rotation angle as a right-handed rotation around + the axis defined by ``axis_tilt`` and ``axis_azimuth``. For example, + with ``axis_tilt=0`` and ``axis_azimuth=180``, ``tracker_theta > 0`` + results in ``surface_azimuth`` to the West while ``tracker_theta < 0`` + results in ``surface_azimuth`` to the East. [degree] + phi : numeric + Projected solar zenith angle. [degrees] + gcr : numeric + Ground coverage ratio, which is the ratio of row slant length to row + spacing (pitch). [unitless] + x0 : numeric, default 0. + Position on the row's slant length, as a fraction of the slant length. + ``x0=0`` corresponds to the bottom of the row. ``x0`` should be less + than ``x1``. [unitless] + x1 : numeric, default 1. + Position on the row's slant length, as a fraction of the slant length. + ``x1`` should be greater than ``x0``. [unitless] + + Returns + ------- + f_x : numeric + Fraction of row slant height from the bottom that is shaded from + direct irradiance. + + References + ---------- + .. [1] Kevin Anderson and Mark Mikofski, "Slope-Aware Backtracking for + Single-Axis Trackers", Technical Report NREL/TP-5K00-76626, July 2020. + https://www.nrel.gov/docs/fy20osti/76626.pdf + """ + # keep track of scalar inputs so that we can have output match at the end + squeeze = [] + if np.isscalar(x0) and np.isscalar(x1): + squeeze.append(0) + if np.isscalar(tracker_rotation): + squeeze.append(1) + + # note: ground slope is already accounted for in phi and gcr, so don't + # apply it here. + # also, we have PSZA instead of solar position, so use fake azimuths to + # trick shaded_fraction1d into accepting it as-is. + # direction of positive phi by right-hand rule: + f_s = shaded_fraction1d(phi, solar_azimuth=90, + axis_azimuth=0, + shaded_row_rotation=tracker_rotation, + collector_width=1, pitch=1/gcr) + + # dimensions: row segment, time + f_s = np.atleast_1d(f_s)[np.newaxis, :] + x0 = np.atleast_1d(x0)[:, np.newaxis] + x1 = np.atleast_1d(x1)[:, np.newaxis] + + swap = tracker_rotation < 0 + x0, x1 = np.where(swap, 1 - x1, x0), np.where(swap, 1 - x0, x1) + + f_s = np.clip((f_s - x0) / (x1 - x0), a_min=0, a_max=1) + f_s = f_s.squeeze(axis=tuple(squeeze)) + + return f_s + + +def _ants2d_singleside(tracker_rotation, cos_aoi, phi, gcr, height, pitch, + dni, dhi, ground_irradiance, albedo, x0, x1, g0, g1, + max_rows): + r""" + Calculate plane-of-array irradiance components on one side of a row + of modules. + + Parameters + ---------- + tracker_rotation : numeric + Tracker rotation angle as a right-handed rotation around + the axis defined by ``axis_tilt`` and ``axis_azimuth``. For example, + with ``axis_tilt=0`` and ``axis_azimuth=180``, ``tracker_theta > 0`` + results in ``surface_azimuth`` to the West while ``tracker_theta < 0`` + results in ``surface_azimuth`` to the East. [degree] + cos_aoi : numeric + Cosine of the angle of incidence of beam irradiance; can be + calculated using :py:func:`pvlib.irradiance.aoi_projection`. [unitless] + phi : numeric + Project solar zenith angle; calculate with + :py:func:`pvlib.shading.projected_solar_zenith_angle`. [degree] + gcr : float + Ground coverage ratio, ratio of row slant length to row spacing. + [unitless] + height : float + Height of the center point of the row above the ground; must be in the + same units as ``pitch``. + pitch : float + Distance between two rows; must be in the same units as ``height``. + dni : numeric + Direct normal irradiance. [Wm⁻²] + dhi : numeric + Diffuse horizontal irradiance. [Wm⁻²] + ground_irradiance : numeric + Irradiance incident on the ground surface, partitioned according + to ``x0`` and ``x1``. Sum of direct and diffuse components. [Wm⁻²] + albedo : numeric + Surface albedo. If a scalar, it is applied to all ground segments and + timestamps. Otherwise, must be specified as an array with shape + (n_ground_segments, n_timestamps). [unitless] + x0 : numeric, default 0 + Position on the row's slant length, as a fraction of the slant length. + ``x0=0`` corresponds to the left side of the row. + ``x0`` should be less than ``x1``. If specified as array, it + must have the same length as ``x1``. [unitless] + x1 : numeric, default 1 + Position on the row's slant length, as a fraction of the slant length. + ``x1=1`` corresponds to the right side of the row. + ``x1`` should be greater than ``x0``. If specified as array, it + must have the same length as ``x0``. [unitless] + g0 : numeric + Position on the ground surface, as a fraction of the row-to-row + spacing. ``g0=0`` corresponds to ground underneath the middle of the + left row. ``g0`` should be less than ``g1``. If specified as array, it + must have the same length as ``g1``. [unitless] + g1 : numeric + Position on the ground surface, as a fraction of the row-to-row + spacing. ``g1=1`` corresponds to ground underneath the middle of the + right row. ``g1`` should be greater than ``g0``. If specified as + array, it must have the same length as ``g0``. [unitless] + max_rows : int + Number of array units (sky wedges, ground segments, etc) to consider. + [unitless] + + Returns + ------- + output : dict of ``numpy.ndarray`` + + ``output`` includes the following quantities: + + - ``poa_global``: total POA irradiance. [Wm⁻²] + - ``poa_diffuse``: total diffuse POA irradiance from all sources. + [Wm⁻²] + - ``poa_direct``: direct POA irradiance. [Wm⁻²] + - ``poa_sky_diffuse``: sky diffuse POA irradiance. [Wm⁻²] + - ``poa_ground_diffuse``: ground-reflected diffuse POA irradiance. + [Wm⁻²] + - ``shaded_fraction``: fraction of row slant height from the bottom + that is shaded from direct irradiance by adjacent rows. [unitless] + + Each array has shape (len(x0), len(tracker_rotation)). + + References + ---------- + .. [1] K. S. Anderson, A. R. Jensen, and C. W. Hansen, "A Bifacial View + Factor Model Considering Terrain Slope and Nonuniform Albedo," + IEEE JPV, 2026. :doi:`10.1109/JPHOTOV.2026.3677506` + """ + # reminder of base dimensions: ground segment, row segment, time + + # in-plane beam component + projection = np.array(np.clip(cos_aoi, a_min=0, a_max=None)) + row_shaded_fraction = _shaded_fraction(tracker_rotation, phi, gcr, x0, x1) + poa_direct = dni * projection * (1 - row_shaded_fraction) + poa_direct = poa_direct[0] # drop ground segment dimension + + # in-plane sky diffuse component + vf_row_sky = utils.vf_row_sky_2d_integ(tracker_rotation, gcr, x0, x1) + poa_sky_diffuse = vf_row_sky * dhi + poa_sky_diffuse = poa_sky_diffuse[0] # drop ground segment dimension + + # in-plane ground-reflected component + vf_row_ground = utils.vf_row_ground_2d_integ(surface_tilt=tracker_rotation, + gcr=gcr, height=height, + pitch=pitch, + x0=x0, x1=x1, g0=g0, g1=g1, + max_rows=max_rows) + poa_ground_diffuse = vf_row_ground * albedo * ground_irradiance + # sum over ground segments + poa_ground_diffuse = np.sum(poa_ground_diffuse, axis=0) + + # add sky and ground-reflected irradiance on the row by irradiance + # component + poa_diffuse = poa_ground_diffuse + poa_sky_diffuse + poa_global = poa_direct + poa_diffuse + + # all arrays are now 2D with shape (n_row_segments, len(tracker_rotation)) + output = { + 'poa_global': poa_global, + 'poa_direct': poa_direct, + 'poa_diffuse': poa_diffuse, + 'poa_sky_diffuse': poa_sky_diffuse, + 'poa_ground_diffuse': poa_ground_diffuse, + 'shaded_fraction': row_shaded_fraction + } + return output + + +def _apply_sky_diffuse_model(dni, dhi, model, solar_zenith, solar_azimuth, + dni_extra, airmass): + + if model in ['haydavies', 'perez']: + # determine circumsolar irradiance, add it to DNI + + if model == 'haydavies': + if dni_extra is None: + raise ValueError(f'Must supply dni_extra for {model} model') + diffuse_model_func = haydavies + extra_kwargs = {} + + elif model == 'perez': + # note: horizon brightening is ignored + if dni_extra is None or airmass is None: + raise ValueError( + f'Must supply dni_extra and airmass for {model} model') + diffuse_model_func = perez + extra_kwargs = {'airmass': airmass} + + kwargs = dict( + dhi=dhi, dni=dni, dni_extra=dni_extra, + solar_zenith=solar_zenith, solar_azimuth=solar_azimuth, + return_components=True + ) + # Call the model first time within the horizontal plane - to subtract + # circumsolar_horizontal from DHI + sky_diffuse_comps_horizontal = diffuse_model_func( + surface_tilt=0, surface_azimuth=180, **kwargs, **extra_kwargs) + circumsolar_horizontal = sky_diffuse_comps_horizontal['poa_circumsolar'] + + # Call the model a second time where circumsolar_normal is facing + # directly towards sun, and can be added to DNI + sky_diffuse_comps_normal = diffuse_model_func( + surface_tilt=solar_zenith, surface_azimuth=solar_azimuth, + **kwargs, **extra_kwargs) + circumsolar_normal = sky_diffuse_comps_normal['poa_circumsolar'] + + dhi = dhi - circumsolar_horizontal + dni = dni + circumsolar_normal + elif model != 'isotropic': + raise ValueError(f"Invalid model: {model}") + + return dni, dhi + + +def _apply_ground_slope(height, pitch, gcr, tracker_rotation, ghi, dni, dhi, + solar_zenith, solar_azimuth, axis_tilt, axis_azimuth, + cross_axis_slope): + slope_azimuth = axis_azimuth + np.degrees( + np.arctan2(sind(cross_axis_slope), + cosd(cross_axis_slope) * sind(axis_tilt)) + ) + slope_tilt = acosd(cosd(axis_tilt) * cosd(cross_axis_slope)) + + height = height * cosd(slope_tilt) + pitch = pitch / cosd(cross_axis_slope) + gcr = gcr * cosd(cross_axis_slope) + tracker_rotation = tracker_rotation - cross_axis_slope + # put back to [-180, 180]: + tracker_rotation = ((tracker_rotation + 180) % 360) - 180 + + ghi = dhi + dni * np.maximum( + aoi_projection(slope_tilt, slope_azimuth, + solar_zenith, solar_azimuth), + 0) + # dhi: no need to adjust; the blocked view is only near the + # the horizon, and that part of the sky is blocked by rows anyway + # dni: no adjustment needed; the measurement plane is not affected + return height, pitch, gcr, tracker_rotation, ghi + + +def get_irradiance(tracker_rotation, axis_azimuth, solar_zenith, solar_azimuth, + gcr, height, pitch, ghi, dhi, dni, + albedo, model='perez', dni_extra=None, airmass=None, + row_segments=1, ground_segments=10, axis_tilt=0, + cross_axis_slope=0, max_rows=None, + return_ground_components=False): + """ + Get front and rear irradiance using the ANTS-2D bifacial irradiance model. + + The ANTS-2D model [1] assumes the PV system comprises parallel, + evenly spaced rows on flat or uniformly sloped ground. Rows can be on fixed + racking or single axis trackers. The model calculates irradiance at a + location far from the ends of any rows, in effect, assuming that the + rows (sheds) are infinitely long. + + The model accounts for the following effects: + + - restricted view of the sky from module surfaces due to the nearby rows. + - restricted view of the ground from module surfaces due to nearby rows. + - restricted view of the sky from the ground due to rows. + - shading of module surfaces by nearby rows. + - nonuniform ground albedo. + - sloped ground surface. + + Parameters + ---------- + tracker_rotation : numeric + Tracker rotation angle as a right-handed rotation around + the axis defined by ``axis_tilt`` and ``axis_azimuth``. For example, + with ``axis_tilt=0`` and ``axis_azimuth=180``, ``tracker_theta > 0`` + results in ``surface_azimuth`` to the West while ``tracker_theta < 0`` + results in ``surface_azimuth`` to the East. [degree] + axis_azimuth : numeric + Axis azimuth angle in degrees. + North = 0°; East = 90°; South = 180°; West = 270° + solar_zenith : numeric + Refraction-corrected solar zenith angle. [degree] + solar_azimuth : numeric + Solar azimuth angle. [degree] + gcr : float + Ground coverage ratio, ratio of row slant length to row spacing. + [unitless] + height : float + Height of the center point of the row above the ground; must be in the + same units as ``pitch``. + pitch : float + Distance between two rows; must be in the same units as ``height``. + ghi : numeric + Global horizontal irradiance. [Wm⁻²] + dhi : numeric + Diffuse horizontal irradiance. [Wm⁻²] + dni : numeric + Direct normal irradiance. [Wm⁻²] + albedo : numeric + Surface albedo. If a scalar, it is applied to all ground segments and + timestamps. Otherwise, must be specified as an array with shape + (``n_ground_segments``, ``len(tracker_rotation)``). [unitless] + model : str, default 'perez' + Irradiance model - can be one of 'isotropic', 'haydavies', or 'perez'. + dni_extra : numeric, optional + Extraterrestrial direct normal irradiance. Required when + ``model='haydavies'`` or ``model='perez'``. [Wm⁻²] + airmass : numeric, optional + Relative airmass. Required when ``model='perez'``. [unitless] + row_segments : int or list of pairs, default 1 + If ``row_segments`` is an int, it defines the number of equal-length + segments the row width is divided into. Otherwise, it must be a list + of pairs ``(x0, x1)`` where ``x0`` and ``x1`` are fractions of the + row width and ``x0 < x1``. Irradiance will be computed and returned + for each segment. + ground_segments : int or list of pairs, default 10 + If ``ground_segments`` is an int, it defines the number of equal-length + segments the ground surface is divided into. Otherwise, it must be + a list of pairs ``(g0, g1)`` where ``g0`` and ``g1`` are fractions of + the pitch and ``g0 < g1``. The pairs must be non-overlapping and must + cover the entire ground surface. ``albedo`` can be specified for + each segment. + axis_tilt : numeric, default 0 + Tilt of the axis of rotation with respect to horizontal. [degree] + cross_axis_slope : numeric, default 0 + The angle, relative to horizontal, of the line formed by the + intersection between the slope containing the tracker axes and a plane + perpendicular to the tracker axes. The cross-axis slope should be + specified using a right-handed convention. For example, trackers with + axis azimuth of 180 degrees (heading south) will have a negative + cross-axis tilt if the tracker axes plane slopes down to the east and + positive cross-axis slope if the tracker axes plane slopes down to the + west. Use :func:`~pvlib.tracking.calc_cross_axis_tilt` to calculate + ``cross_axis_slope``. [degrees] + max_rows : int, optional + Number of array units (sky wedges, ground segments, etc) to consider. + If not specified, units will be considered to within 4 degrees of the + horizon. [unitless] + return_ground_components : bool, default False + If True, also return the direct and diffuse irradiance incident on the + ground. These values are returned in a second dict. + + Returns + ------- + output : dict or DataFrame + ``output`` is a DataFrame when inputs are Series + and ``row_segments=1``, a dict of scalars when inputs are scalars + and ``row_segments=1``, and a dict of ``np.ndarray`` + otherwise. The following quantities are included: + + - ``poa_global``: sum of front- and back-side incident irradiance. + [Wm⁻²] + - ``poa_front``: total incident irradiance on the front surface. [Wm⁻²] + - ``poa_back``: total incident irradiance on the back surface. [Wm⁻²] + - ``poa_front_direct``: direct irradiance incident on the front + surface. [Wm⁻²] + - ``poa_front_diffuse``: total diffuse irradiance incident on the front + surface. [Wm⁻²] + - ``poa_front_sky_diffuse``: sky diffuse irradiance incident on the + front surface. [Wm⁻²] + - ``poa_front_ground_diffuse``: ground-reflected diffuse irradiance + incident on the front surface. [Wm⁻²] + - ``shaded_fraction_front``: fraction of row slant height that is + shaded from direct irradiance on the front surface by adjacent + rows. [unitless] + - ``poa_back_direct``: direct irradiance incident on the back + surface. [Wm⁻²] + - ``poa_back_diffuse``: total diffuse irradiance incident on the back + surface. [Wm⁻²] + - ``poa_back_sky_diffuse``: sky diffuse irradiance incident on the + back surface. [Wm⁻²] + - ``poa_back_ground_diffuse``: ground-reflected diffuse irradiance + incident on the back surface. [Wm⁻²] + - ``shaded_fraction_back``: fraction of row slant height that is + shaded from direct irradiance on the back surface by adjacent + rows. [unitless] + + ground_irradiance : dict or DataFrame + ``ground_irradiance`` is a DataFrame when inputs are Series + and ``n_ground_segments=1``, a dict of scalars when inputs are scalars + and ``n_ground_segments=1``, and a dict of ``np.ndarray`` + otherwise. Only returned when ``return_ground_components=True``. + The following quantities are included: + + - ``ground_direct``: direct irradiance incident on the ground surface. + [Wm⁻²] + - ``ground_diffuse``: diffuse irradiance incident on the ground + surface. [Wm⁻²] + + References + ---------- + .. [1] K. S. Anderson, A. R. Jensen, and C. W. Hansen, "A Bifacial View + Factor Model Considering Terrain Slope and Nonuniform Albedo," + IEEE JPV, 2026. :doi:`10.1109/JPHOTOV.2026.3677506` + """ + + # so we can return scalars out if needed + maybe_array_inputs = [ + tracker_rotation, axis_azimuth, solar_zenith, solar_azimuth, + gcr, height, pitch, ghi, dhi, dni, + albedo, dni_extra, airmass, axis_tilt, cross_axis_slope] + all_scalar_inputs = all([ + np.isscalar(x) or x is None for x in maybe_array_inputs + ]) + try: + # get the index of the first pandas input, if there is one + pd_index = next( + x.index for x in maybe_array_inputs if isinstance(x, pd.Series) + ) + except StopIteration: + pd_index = None # no pandas inputs + + # preparation steps + + dni, dhi = _apply_sky_diffuse_model(dni, dhi, model, solar_zenith, + solar_azimuth, dni_extra, airmass) + + true_tracker_rotation = tracker_rotation + + if axis_tilt != 0 or cross_axis_slope != 0: + height, pitch, gcr, tracker_rotation, ghi = _apply_ground_slope( + height, pitch, gcr, tracker_rotation, ghi, dni, dhi, + solar_zenith, solar_azimuth, axis_tilt, axis_azimuth, + cross_axis_slope + ) + + if np.isscalar(row_segments): + x_row = np.linspace(0, 1, row_segments+1) + x0, x1 = x_row[:-1], x_row[1:] + else: + x0 = np.array([pair[0] for pair in row_segments]) + x1 = np.array([pair[1] for pair in row_segments]) + + if np.isscalar(ground_segments): + x_ground = np.linspace(0, 1, ground_segments+1) + g0, g1 = x_ground[:-1], x_ground[1:] + else: + g0 = np.array([pair[0] for pair in ground_segments]) + g1 = np.array([pair[1] for pair in ground_segments]) + + # dimensions: ground segment, row segment, time + albedo = np.atleast_2d(albedo)[:, np.newaxis, :] + ghi = np.atleast_1d(ghi)[np.newaxis, np.newaxis, :] + dhi = np.atleast_1d(dhi)[np.newaxis, np.newaxis, :] + dni = np.atleast_1d(dni)[np.newaxis, np.newaxis, :] + tracker_rotation = np.atleast_1d(tracker_rotation) + + # Calculate some geometric quantities + # rows to consider in front and behind current row + # ensures that view factors to the sky are computed to within 4 degrees + # of the horizon + if max_rows is None: + max_rows = np.ceil(height / (pitch * tand(4))) + + phi = projected_solar_zenith_angle(solar_zenith, solar_azimuth, + axis_tilt, axis_azimuth) + phi = phi - cross_axis_slope + + # compute this here, as it is expensive and does not differ between the + # front and back sides + vf_gnd_sky = utils.vf_ground_sky_2d_integ( + tracker_rotation, gcr, height, pitch, g0=g0, g1=g1, max_rows=max_rows) + vf_gnd_sky = vf_gnd_sky[:, np.newaxis, :] + + # irradiance components incident on ground surface + ground_unshaded_fraction = utils._unshaded_ground_fraction( + tracker_rotation, phi, gcr, + pitch=pitch, height=height, g0=g0, g1=g1, max_rows=max_rows) + + ground_shaded_fraction = 1 - ground_unshaded_fraction + ground_shaded_fraction = ground_shaded_fraction[:, np.newaxis, :] + + ground_direct = (1-ground_shaded_fraction) * (ghi - dhi) + ground_diffuse = vf_gnd_sky * dhi + ground_total = ground_direct + ground_diffuse + + # inputs shared between front and back calculations + params = dict(phi=phi, gcr=gcr, height=height, pitch=pitch, dni=dni, + dhi=dhi, ground_irradiance=ground_total, albedo=albedo, + g0=g0, g1=g1, max_rows=max_rows) + + # front + front_orientation = calc_surface_orientation(true_tracker_rotation, + axis_tilt, axis_azimuth) + cos_aoi_front = aoi_projection(**front_orientation, + solar_zenith=solar_zenith, + solar_azimuth=solar_azimuth) + poa_front = _ants2d_singleside(tracker_rotation, cos_aoi_front, + x0=x0, x1=x1, **params) + + # back + tracker_rotation_back = true_tracker_rotation + 180 + tracker_rotation_back = ((tracker_rotation_back + 180) % 360) - 180 + back_orientation = calc_surface_orientation(tracker_rotation_back, + axis_tilt, axis_azimuth) + cos_aoi_back = aoi_projection(**back_orientation, + solar_zenith=solar_zenith, + solar_azimuth=solar_azimuth) + tracker_rotation_back = tracker_rotation + 180 + tracker_rotation_back = ((tracker_rotation_back + 180) % 360) - 180 + poa_back = _ants2d_singleside(tracker_rotation_back, cos_aoi_back, + x0=1-x1, x1=1-x0, **params) + + colmap_front = { + 'poa_global': 'poa_front', + 'poa_direct': 'poa_front_direct', + 'poa_diffuse': 'poa_front_diffuse', + 'poa_sky_diffuse': 'poa_front_sky_diffuse', + 'poa_ground_diffuse': 'poa_front_ground_diffuse', + 'shaded_fraction': 'shaded_fraction_front', + } + colmap_back = { + k: v.replace("front", "back") for k, v in colmap_front.items() + } + for old_key, new_key in colmap_front.items(): + poa_front[new_key] = poa_front.pop(old_key) + for old_key, new_key in colmap_back.items(): + poa_back[new_key] = poa_back.pop(old_key) + out = {**poa_front, **poa_back} + + if row_segments == 1: + for k, v in out.items(): + out[k] = v[0] # drop row segment dimension + + if all_scalar_inputs: + # drop the second dimension too, so scalars are returned + for k, v in out.items(): + out[k] = float(v[0]) + + elif pd_index is not None: + out = pd.DataFrame(out, index=pd_index) + + if return_ground_components: + squeeze = [] + if ground_segments == 1: + squeeze.append(0) # drop ground segment dimension + squeeze.append(1) # always drop the row segment dimension + if all_scalar_inputs: + squeeze.append(2) # drop time dimension + squeeze = tuple(squeeze) + out_ground = { + 'ground_direct': ground_direct.squeeze(axis=squeeze), + 'ground_diffuse': ground_diffuse.squeeze(axis=squeeze), + } + if squeeze == (0, 1, 2): + for k, v in out_ground.items(): + out_ground[k] = float(v) + + if pd_index is not None and squeeze == (0, 1): + out_ground = pd.DataFrame(out_ground, index=pd_index) + + return out, out_ground + + return out diff --git a/pvlib/bifacial/infinite_sheds.py b/pvlib/bifacial/infinite_sheds.py index f99c536315..15ef4a89bf 100644 --- a/pvlib/bifacial/infinite_sheds.py +++ b/pvlib/bifacial/infinite_sheds.py @@ -4,7 +4,7 @@ import numpy as np import pandas as pd -from pvlib.tools import cosd, sind, tand +from pvlib.tools import cosd, sind, tand, atand from pvlib.bifacial import utils from pvlib.irradiance import beam_component, aoi, haydavies @@ -90,7 +90,7 @@ def _poa_sky_diffuse_pv(dhi, gcr, surface_tilt): poa_sky_diffuse_pv : numeric Total sky diffuse irradiance incident on the PV surface. [W/m^2] """ - vf_integ = utils.vf_row_sky_2d_integ(surface_tilt, gcr, 0., 1.) + vf_integ = utils.vf_row_sky_2d_integ(surface_tilt, gcr, x0=0., x1=1.) return dhi * vf_integ @@ -117,7 +117,8 @@ def _poa_ground_pv(poa_ground, gcr, surface_tilt): numeric Ground diffuse irradiance on the row plane. [W/m^2] """ - vf_integ = utils.vf_row_ground_2d_integ(surface_tilt, gcr, 0., 1.) + vf_integ = utils.vf_row_ground_2d_integ(surface_tilt, gcr, x0=0., x1=1., + g0=0., g1=1.) return poa_ground * vf_integ @@ -322,15 +323,18 @@ def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith, max_rows = np.ceil(height / (pitch * tand(5))) # fraction of ground between rows that is illuminated accounting for # shade from panels. [1], Eq. 4 - f_gnd_beam = utils._unshaded_ground_fraction( - surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, gcr) + tanphi = utils._solar_projection_tangent(solar_zenith, solar_azimuth, + surface_azimuth) + phi = atand(tanphi) + f_gnd_beam = utils._unshaded_ground_fraction(surface_tilt, phi, gcr) + # integrated view factor from the ground to the sky, integrated between # adjacent rows interior to the array # method differs from [1], Eq. 7 and Eq. 8; height is defined at row # center rather than at row lower edge as in [1]. vf_gnd_sky = utils.vf_ground_sky_2d_integ( - surface_tilt, gcr, height, pitch, max_rows, npoints, - vectorize) + surface_tilt, gcr, height, pitch, g0=0, g1=1, max_rows=max_rows, + npoints=npoints, vectorize=vectorize) # fraction of row slant height that is shaded from direct irradiance f_x = _shaded_fraction(solar_zenith, solar_azimuth, surface_tilt, surface_azimuth, gcr) diff --git a/pvlib/bifacial/utils.py b/pvlib/bifacial/utils.py index a69a3e18d4..eb6a5ec7de 100644 --- a/pvlib/bifacial/utils.py +++ b/pvlib/bifacial/utils.py @@ -4,7 +4,8 @@ """ import numpy as np from pvlib.tools import sind, cosd, tand -from scipy.integrate import trapezoid +import warnings +from pvlib._deprecation import pvlibDeprecationWarning def _solar_projection_tangent(solar_zenith, solar_azimuth, surface_azimuth): @@ -37,36 +38,42 @@ def _solar_projection_tangent(solar_zenith, solar_azimuth, surface_azimuth): return tan_phi -def _unshaded_ground_fraction(surface_tilt, surface_azimuth, solar_zenith, - solar_azimuth, gcr, max_zenith=87): +def _unshaded_ground_fraction(surface_tilt, phi, gcr, height=None, + pitch=None, g0=0, g1=1, max_rows=10, + max_zenith=85): r""" Calculate the fraction of the ground with incident direct irradiance. - .. math:: - F_{gnd,sky} = 1 - \min{\left(1, \text{GCR} \left|\cos \beta + - \sin \beta \tan \phi \right|\right)} - - where :math:`\beta` is the surface tilt and :math:`\phi` is the angle - from vertical of the sun vector projected to a vertical plane that - contains the row azimuth `surface_azimuth`. - Parameters ---------- surface_tilt : numeric Surface tilt angle. The tilt angle is defined as degrees from horizontal, e.g., surface facing up = 0, surface facing horizon = 90. [degree] - surface_azimuth : numeric - Azimuth of the module surface, i.e., North=0, East=90, South=180, - West=270. [degree] - solar_zenith : numeric - Solar zenith angle. [degree]. - solar_azimuth : numeric - Solar azimuth. [degree]. + phi : numeric + Projected solar zenith angle. [degree]. gcr : float Ground coverage ratio, which is the ratio of row slant length to row spacing (pitch). [unitless] - max_zenith : numeric, default 87 + height : float, optional + Height of the center point of the row above the ground; must be in the + same units as ``pitch``. Required if ``g0`` is not zero or ``g1`` is + not one. + pitch : float, optional + Distance between two rows; must be in the same units as ``height``. + Required if ``g0`` is not zero or ``g1`` is not one. + g0 : numeric, default 0 + Position on the ground surface, as a fraction of the row-to-row + spacing. ``g0=0`` corresponds to ground underneath the middle of the + left row. ``g0`` should be less than ``g1``. [unitless] + g1 : numeric, default 1 + Position on the ground surface, as a fraction of the row-to-row + spacing. ``g1=1`` corresponds to ground underneath the middle of the + right row. ``g1`` should be greater than ``g0``. [unitless] + max_rows : int, default 10 + Maximum number of rows to consider on either side of the current + row. [unitless] + max_zenith : numeric, default 85 Maximum zenith angle. For solar_zenith > max_zenith, unshaded ground fraction is set to 0. [degree] @@ -83,13 +90,70 @@ def _unshaded_ground_fraction(surface_tilt, surface_azimuth, solar_zenith, Photovoltaic Specialists Conference (PVSC), 2019, pp. 1282-1287. :doi:`10.1109/PVSC40753.2019.8980572`. """ - tan_phi = _solar_projection_tangent(solar_zenith, solar_azimuth, - surface_azimuth) - f_gnd_beam = 1.0 - np.minimum( - 1.0, gcr * np.abs(cosd(surface_tilt) + sind(surface_tilt) * tan_phi)) - # [1], Eq. 4 - f_gnd_beam = np.where(solar_zenith > max_zenith, 0., f_gnd_beam) - return f_gnd_beam # 1 - min(1, abs()) < 1 always + if np.isscalar(g0) and g0 == 0 and np.isscalar(g1) and g1 == 1: + # height and pitch have no effect, so set to arbitrary values so + # that they can be optional parameters + height = 1 + pitch = 1 / gcr + + # keep track of scalar inputs so that we can have output match at the end + squeeze = [] + if np.isscalar(g0) and np.isscalar(g1): + squeeze.append(0) + if np.isscalar(surface_tilt): + squeeze.append(1) + + # dimensions: k/max_rows, ground segment, time + + surface_tilt = np.atleast_1d(surface_tilt)[np.newaxis, np.newaxis, :] + phi = np.atleast_1d(phi)[np.newaxis, np.newaxis, :] + + g0 = np.atleast_1d(g0)[np.newaxis, :, np.newaxis] + g1 = np.atleast_1d(g1)[np.newaxis, :, np.newaxis] + + # TODO seems like this should be np.arange(-max_rows, max_rows+1)? + # see GH #1867 + k = np.arange(-max_rows, max_rows)[:, np.newaxis, np.newaxis] + + collector_width = pitch * gcr + Lcostheta = collector_width * cosd(surface_tilt) + Lsintheta = collector_width * sind(surface_tilt) + tanphi = tand(phi) + + # a, b: boundaries of ground segment + # d, c: left/right shading module edges + c = (k*pitch + 0.5 * Lcostheta, height + 0.5 * Lsintheta) + d = (k*pitch - 0.5 * Lcostheta, height - 0.5 * Lsintheta) + + cp = c[0] + c[1] * tanphi + dp = d[0] + d[1] * tanphi + swap = dp > cp + cp, dp = np.where(swap, dp, cp), np.where(swap, cp, dp) + + a = g0*pitch + b = g1*pitch + + # individual contributions from all k rows + fs = np.full_like(cp, 1.0) + # fs = np.where((dp <= a) & (cp >= b), 1.0, fs) # fill value already 1.0 + fs = np.where((dp <= a) & (a <= cp) & (cp < b), (cp - a) / (b - a), fs) + fs = np.where((dp < a) & (cp < a), 0.0, fs) + fs = np.where((a < dp) & (dp <= b) & (cp >= b), (b - dp) / (b - a), fs) + fs = np.where((a < dp) & (dp < b) & (a < cp) & (cp < b), + (cp - dp) / (b - a), fs) + fs = np.where((dp > b) & (cp > b), 0.0, fs) + + # total shaded fraction is sum of individuals; note that shadows + # never overlap in this model, except when shaded fraction is 100% anyway + f_gnd_beam = 1 - np.clip(np.sum(fs, axis=0), 0, 1) # sum along k dimension + + # using phi is more convenient, and I think better, than using zenith + phi = phi[0, :, :] # drop k dimension for the next line + f_gnd_beam = np.where(np.abs(phi) > max_zenith, 0., f_gnd_beam) + + # dimensions are now ground_segment, time + f_gnd_beam = f_gnd_beam.squeeze(axis=tuple(squeeze)) + return f_gnd_beam def vf_ground_sky_2d(rotation, gcr, x, pitch, height, max_rows=10): @@ -111,11 +175,11 @@ def vf_ground_sky_2d(rotation, gcr, x, pitch, height, max_rows=10): Position on the ground between two rows, as a fraction of the pitch. x = 0 corresponds to the point on the ground directly below the center point of a row. Positive x is towards the right. [unitless] + pitch : float + Distance between two rows; must be in the same units as ``height``. height : float Height of the center point of the row above the ground; must be in the same units as ``pitch``. - pitch : float - Distance between two rows; must be in the same units as ``height``. max_rows : int, default 10 Maximum number of rows to consider on either side of the current row. [unitless] @@ -174,8 +238,8 @@ def vf_ground_sky_2d(rotation, gcr, x, pitch, height, max_rows=10): return vf -def vf_ground_sky_2d_integ(surface_tilt, gcr, height, pitch, max_rows=10, - npoints=100, vectorize=False): +def vf_ground_sky_2d_integ(tracker_rotation, gcr, height, pitch, g0=0, g1=1, + max_rows=10, npoints=None, vectorize=None): """ Integrated view factor to the sky from the ground underneath interior rows of the array. @@ -192,13 +256,29 @@ def vf_ground_sky_2d_integ(surface_tilt, gcr, height, pitch, max_rows=10, same units as ``pitch``. pitch : float Distance between two rows. Must be in the same units as ``height``. + g0 : numeric, default 0 + Position on the ground surface, as a fraction of the row-to-row + spacing. ``g0=0`` corresponds to ground underneath the middle of the + left row. ``g0`` should be less than ``g1``. [unitless] + g1 : numeric, default 1 + Position on the ground surface, as a fraction of the row-to-row + spacing. ``g1=1`` corresponds to ground underneath the middle of the + right row. ``g1`` should be greater than ``g0``. [unitless] max_rows : int, default 10 Maximum number of rows to consider in front and behind the current row. - npoints : int, default 100 - Number of points used to discretize distance along the ground. - vectorize : bool, default False - If True, vectorize the view factor calculation across ``surface_tilt``. - This increases speed with the cost of increased memory usage. + npoints : int, optional + + .. deprecated:: TODO + + This parameter has no effect; integrated view factors are now + calculated exactly instead of with discretized approximations. + + vectorize : bool, optional + + .. deprecated:: TODO + + This parameter has no effect; calculations are now vectorized + with no memory usage penality. Returns ------- @@ -206,23 +286,78 @@ def vf_ground_sky_2d_integ(surface_tilt, gcr, height, pitch, max_rows=10, Integration of view factor over the length between adjacent, interior rows. Shape matches that of ``surface_tilt``. [unitless] """ - # Abuse vf_ground_sky_2d by supplying surface_tilt in place - # of a signed rotation. This is OK because - # 1) z span the full distance between 2 rows, and - # 2) max_rows is set to be large upstream, and - # 3) _vf_ground_sky_2d considers [-max_rows, +max_rows] - # The VFs to the sky will thus be symmetric around z=0.5 - z = np.linspace(0, 1, npoints) - rotation = np.atleast_1d(surface_tilt) - if vectorize: - fz_sky = vf_ground_sky_2d(rotation, gcr, z, pitch, height, max_rows) - else: - fz_sky = np.zeros((npoints, len(rotation))) - for k, r in enumerate(rotation): - vf = vf_ground_sky_2d(r, gcr, z, pitch, height, max_rows) - fz_sky[:, k] = vf[:, 0] # remove spurious rotation dimension - # calculate the integrated view factor for all of the ground between rows - return trapezoid(fz_sky, z, axis=0) + if npoints is not None or vectorize is not None: + msg = ( + "The `npoints` and `vectorize` parameters have no effect and will " + "be removed in a future version." # TODO make this better + ) + warnings.warn(msg, pvlibDeprecationWarning) + + # keep track of scalar inputs so that we can have output match at the end + squeeze = [] + if np.isscalar(g0) and np.isscalar(g1): + squeeze.append(0) + if np.isscalar(tracker_rotation): + squeeze.append(1) + + # dimensions: k/max_rows, ground segment, time + + tracker_rotation = \ + np.atleast_1d(tracker_rotation)[np.newaxis, np.newaxis, :] + + g0 = np.atleast_1d(g0)[np.newaxis, :, np.newaxis] + g1 = np.atleast_1d(g1)[np.newaxis, :, np.newaxis] + + # TODO seems like this should be np.arange(-max_rows, max_rows+1)? + # see GH #1867 + k = np.arange(-max_rows, max_rows)[:, np.newaxis, np.newaxis] + + collector_width = pitch * gcr + Lcostheta = collector_width * cosd(tracker_rotation) + Lsintheta = collector_width * sind(tracker_rotation) + + # primary crossed string points: + # a, b: boundaries of ground segment + # c, d: upper module edges + a = (g0*pitch, 0) + b = (g1*pitch, 0) + sign = np.sign(tracker_rotation) + c = ((k+1)*pitch + sign * 0.5 * Lcostheta, height + sign * 0.5 * Lsintheta) + d = (c[0] - pitch, c[1]) + + # view obstruction points (module edges, but need to figure out which ones) + + # first decide whether the left obstruction is the left or right mod edge + left = (k*pitch - 0.5 * Lcostheta, height - 0.5 * Lsintheta) + right = (k*pitch + 0.5 * Lcostheta, height + 0.5 * Lsintheta) + angle_left = _angle(a, left) + angle_right = _angle(a, right) + ob_left = ( + np.where(angle_left > angle_right, right[0], left[0]), + np.where(angle_left > angle_right, right[1], left[1]) + ) + + # now for the right obstruction + left = (left[0] + pitch, left[1]) + right = (right[0] + pitch, right[1]) + angle_left = _angle(b, left) + angle_right = _angle(b, right) + ob_right = ( + np.where(angle_left > angle_right, left[0], right[0]), + np.where(angle_left > angle_right, left[1], right[1]) + ) + + # hottel string lengths, considering obstructions + ac, ad, bc, bd = _obstructed_string_lengths(a, b, c, d, ob_left, ob_right) + + # crossed string formula for VF + vf_slats = 0.5 * (1/((g1 - g0) * pitch)) * ((ac + bd) - (bc + ad)) + vf_total = np.sum(np.maximum(vf_slats, 0), axis=0) # sum along k dimension + + # dimensions are now ground_segment, row_segment, time + vf_total = vf_total.squeeze(axis=tuple(squeeze)) + + return vf_total def _vf_poly(surface_tilt, gcr, x, delta): @@ -298,11 +433,11 @@ def vf_row_sky_2d_integ(surface_tilt, gcr, x0=0, x1=1): Ratio of the row slant length to the row spacing (pitch). [unitless] x0 : numeric, default 0 Position on the row's slant length, as a fraction of the slant length. - x0=0 corresponds to the bottom of the row. x0 should be less than x1. - [unitless] + ``x0=0`` corresponds to the bottom of the row. ``x0`` should be less + than ``x1``. [unitless] x1 : numeric, default 1 Position on the row's slant length, as a fraction of the slant length. - x1 should be greater than x0. [unitless] + ``x1`` should be greater than ``x0``. [unitless] Returns ------- @@ -311,6 +446,23 @@ def vf_row_sky_2d_integ(surface_tilt, gcr, x0=0, x1=1): from x0 to x1. [unitless] ''' + # keep track of scalar inputs so that we can have output match at the end + squeeze = [] + if np.isscalar(x0) and np.isscalar(x1): + squeeze.append(0) + if np.isscalar(surface_tilt): + squeeze.append(1) + + # dimensions: row segment, time + + surface_tilt = np.atleast_1d(surface_tilt)[np.newaxis, :] + + x0 = np.atleast_1d(x0)[:, np.newaxis] + x1 = np.atleast_1d(x1)[:, np.newaxis] + + swap = surface_tilt < 0 + x0, x1 = np.where(swap, 1 - x1, x0), np.where(swap, 1 - x0, x1) + u = np.abs(x1 - x0) p0 = _vf_poly(surface_tilt, gcr, 1 - x0, -1) p1 = _vf_poly(surface_tilt, gcr, 1 - x1, -1) @@ -319,6 +471,7 @@ def vf_row_sky_2d_integ(surface_tilt, gcr, x0=0, x1=1): vf_row_sky_2d(surface_tilt, gcr, x0), 0.5*(1 + 1/u * (p1 - p0)) ) + result = result.squeeze(axis=tuple(squeeze)) return result @@ -351,7 +504,8 @@ def vf_row_ground_2d(surface_tilt, gcr, x): return 0.5 * (1 - (1/gcr * cosd(surface_tilt) + x)/p) -def vf_row_ground_2d_integ(surface_tilt, gcr, x0=0, x1=1): +def vf_row_ground_2d_integ(surface_tilt, gcr, height=None, pitch=None, + x0=0, x1=1, g0=0, g1=1, max_rows=20): r''' Calculate the average view factor to the ground from a segment of the row surface between x0 and x1. @@ -367,13 +521,31 @@ def vf_row_ground_2d_integ(surface_tilt, gcr, x0=0, x1=1): = 0, surface facing horizon = 90. [degree] gcr : numeric Ratio of the row slant length to the row spacing (pitch). [unitless] + height : float, optional + Height of the center point of the row above the ground; must be in the + same units as ``pitch``. Required if ``g0`` or ``x0` is not zero or + if ``g1`` or ``x1`` is not one. + pitch : float, optional + Distance between two rows; must be in the same units as ``height``. + Required if ``g0`` or ``x0` is not zero or if ``g1`` or ``x1`` + is not one. x0 : numeric, default 0. Position on the row's slant length, as a fraction of the slant length. - x0=0 corresponds to the bottom of the row. x0 should be less than x1. - [unitless] + ``x0=0`` corresponds to the bottom of the row. ``x0`` should be less + than ``x1``. [unitless] x1 : numeric, default 1. Position on the row's slant length, as a fraction of the slant length. - x1 should be greater than x0. [unitless] + ``x1`` should be greater than ``x0``. [unitless] + g0 : numeric, default 0 + Position on the ground surface, as a fraction of the row-to-row + spacing. ``g0=0`` corresponds to ground underneath the middle of the + left row. ``g0`` should be less than ``g1``. [unitless] + g1 : numeric, default 1 + Position on the ground surface, as a fraction of the row-to-row + spacing. ``g1=0`` corresponds to ground underneath the middle of the + right row. ``g1`` should be greater than ``g0``. [unitless] + max_rows : int, default 20 + Maximum number of rows to consider in front and behind the current row. Returns ------- @@ -382,12 +554,141 @@ def vf_row_ground_2d_integ(surface_tilt, gcr, x0=0, x1=1): [unitless] ''' - u = np.abs(x1 - x0) - p0 = _vf_poly(surface_tilt, gcr, x0, 1) - p1 = _vf_poly(surface_tilt, gcr, x1, 1) - with np.errstate(divide='ignore'): - result = np.where(u < 1e-6, - vf_row_ground_2d(surface_tilt, gcr, x0), - 0.5*(1 - 1/u * (p1 - p0)) - ) - return result + if all(np.isscalar(x) for x in [x0, x1, g0, g1]) and ( + g0 == 0 and g1 == 1 and x0 == 0 and x1 == 1): + # height and pitch have no effect, so set to arbitrary values so + # that they can be optional parameters + height = 1 + pitch = 1 / gcr + + # keep track of scalar inputs so that we can have output match at the end + squeeze = [] + if np.isscalar(g0) and np.isscalar(g1): + squeeze.append(0) + if np.isscalar(x0) and np.isscalar(x1): + squeeze.append(1) + if np.isscalar(surface_tilt): + squeeze.append(2) + + # dimensions: k/max_rows, ground segment, row segment, time + + # cheat a little to prevent numerical issues with surface_tilt==180, -180 + surface_tilt = np.where(surface_tilt == 180, 179.9999, surface_tilt) + surface_tilt = np.where(surface_tilt == -180, -179.9999, surface_tilt) + + surface_tilt = \ + np.atleast_1d(surface_tilt)[np.newaxis, np.newaxis, np.newaxis, :] + + x0 = np.atleast_1d(x0)[np.newaxis, np.newaxis, :, np.newaxis] + x1 = np.atleast_1d(x1)[np.newaxis, np.newaxis, :, np.newaxis] + g0 = np.atleast_1d(g0)[np.newaxis, :, np.newaxis, np.newaxis] + g1 = np.atleast_1d(g1)[np.newaxis, :, np.newaxis, np.newaxis] + + # TODO seems like this should be np.arange(-max_rows, max_rows+1)? + # see GH #1867 + k = np.arange(-max_rows, max_rows)[:, np.newaxis, np.newaxis, np.newaxis] + + collector_width = pitch * gcr + Lcostheta = collector_width * cosd(surface_tilt) + Lsintheta = collector_width * sind(surface_tilt) + + # view obstruction points (lower module edges) + # use a number slightly larger than 0.5 because the obstruction must + # be a nonzero distance from all points the VF could be calculated from + ob_right = (-pitch - 0.5001 * Lcostheta, height - 0.5001 * abs(Lsintheta)) + ob_left = (ob_right[0] + pitch, ob_right[1]) + + invert = surface_tilt < 0 + temp = ob_right[0] + ob_right = (np.where(invert, -ob_left[0], ob_right[0]), ob_right[1]) + ob_left = (np.where(invert, -temp, ob_left[0]), ob_left[1]) + + # primary crossed string points: + # a, b: positions on module + # c, d: boundaries of ground segment + + a = ((x0-0.5) * Lcostheta, height + (x0-0.5) * Lsintheta) + b = ((x1-0.5) * Lcostheta, height + (x1-0.5) * Lsintheta) + c = ((k+g0)*pitch, 0) + d = ((k+g1)*pitch, 0) + + # hottel string lengths, considering obstructions + ac, ad, bc, bd = _obstructed_string_lengths(a, b, c, d, ob_left, ob_right) + + # crossed string formula for VF + vf_slats = 1 / (2 * (x1 - x0) * collector_width) * ((ac + bd) - (bc + ad)) + vf_total = np.sum(np.maximum(vf_slats, 0), axis=0) # sum along k dimension + + # dimensions are now ground_segment, row_segment, time + vf_total = vf_total.squeeze(axis=tuple(squeeze)) + + return vf_total + + +def _obstructed_string_lengths(a, b, c, d, ob_left, ob_right): + # string length calculations for Hottel's crossed strings method, + # considering view obstructions from the left and right. + # all inputs are (x, y) points. + + # unobstructed (straight-line) distances + dist_ac = _dist(a, c) + dist_ad = _dist(a, d) + dist_bc = _dist(b, c) + dist_bd = _dist(b, d) + dist_a_ob_left = _dist(a, ob_left) + dist_a_ob_right = _dist(a, ob_right) + dist_b_ob_left = _dist(b, ob_left) + dist_b_ob_right = _dist(b, ob_right) + dist_ob_left_c = _dist(ob_left, c) + dist_ob_right_c = _dist(ob_right, c) + dist_ob_left_d = _dist(ob_left, d) + dist_ob_right_d = _dist(ob_right, d) + + # angles + ang_ac = _angle(a, c) + ang_ad = _angle(a, d) + ang_bc = _angle(b, c) + ang_bd = _angle(b, d) + ang_a_ob_left = _angle(a, ob_left) + ang_a_ob_right = _angle(a, ob_right) + ang_b_ob_left = _angle(b, ob_left) + ang_b_ob_right = _angle(b, ob_right) + + # obstructed distances + ac = np.where(ang_ac - ang_a_ob_left > 1e-8, + dist_a_ob_left + dist_ob_left_c, + dist_ac) + ac = np.where((ang_a_ob_right - ang_ac > 1e-8), + dist_a_ob_right + dist_ob_right_c, + ac) + + ad = np.where(ang_ad - ang_a_ob_left > 1e-8, + dist_a_ob_left + dist_ob_left_d, + dist_ad) + ad = np.where(ang_a_ob_right - ang_ad > 1e-8, + dist_a_ob_right + dist_ob_right_d, + ad) + + bc = np.where(ang_bc - ang_b_ob_left > 1e-8, + dist_b_ob_left + dist_ob_left_c, + dist_bc) + bc = np.where(ang_b_ob_right - ang_bc > 1e-8, + dist_b_ob_right + dist_ob_right_c, + bc) + + bd = np.where(ang_bd - ang_b_ob_left > 1e-8, + dist_b_ob_left + dist_ob_left_d, + dist_bd) + bd = np.where(ang_b_ob_right - ang_bd > 1e-8, + dist_b_ob_right + dist_ob_right_d, + bd) + + return ac, ad, bc, bd + + +def _dist(p1, p2): + return ((p1[0] - p2[0])**2 + (p1[1] - p2[1])**2)**0.5 + + +def _angle(p1, p2): + return np.arctan2(p2[1] - p1[1], p2[0] - p1[0]) diff --git a/tests/bifacial/test_ants2d.py b/tests/bifacial/test_ants2d.py new file mode 100644 index 0000000000..280c3e4d6d --- /dev/null +++ b/tests/bifacial/test_ants2d.py @@ -0,0 +1,614 @@ +""" +test ants2d +""" + +import numpy as np +import pandas as pd +import pvlib +from pvlib.bifacial import ants2d + +import pytest +from numpy.testing import assert_allclose + + +def test__shaded_fraction(): + # special angles + tracker_rotation = np.array([60, 60, 60, 60]) + phi = np.array([60, 60, 60, 60]) + gcr = np.array([1, 0.75, 2/3, 0.5]) + expected = np.array([0.5, 1/3, 0.25, 0]) + fs = ants2d._shaded_fraction(tracker_rotation, phi, gcr) + assert_allclose(fs, expected) + fs = ants2d._shaded_fraction(-tracker_rotation, -phi, gcr) + assert_allclose(fs, expected) + + # sun too high for shade + assert 0 == ants2d._shaded_fraction(10, 20, 0.5) + assert 0 == ants2d._shaded_fraction(-10, -20, 0.5) + + # sun behind the modules (AOI > 90) + # (debatable whether this should be zero or one) + assert 0 == ants2d._shaded_fraction(45, -50, 0.5) + assert 0 == ants2d._shaded_fraction(-45, 50, 0.5) + + # edge cases + tracker_rotation = np.array([0, 0, 0, 90, 90, 90, -90, -90, -90]) + phi = np.array([0, 90, -90, 0, 90, -90, 0, 90, -90]) + gcr = 0.5 + # (some of these are debatable as well) + expected = np.array([0, 0, 0, 0, 1, 1, 0, 1, 1]) + fs = ants2d._shaded_fraction(tracker_rotation, phi, gcr) + assert_allclose(fs, expected) + + +def test__shaded_fraction_x0x1(): + fs = ants2d._shaded_fraction(np.array([60, -60]), np.array([60, -60]), + 2/3, x0=[0, 0.5], x1=[0.5, 1]) + assert_allclose(fs, np.array([[0.5, 0.0], [0.0, 0.5]])) + + +@pytest.mark.parametrize('model', ['perez', 'haydavies']) +def test__apply_sky_diffuse_model(model): + inputs = {'dni': 900, 'dhi': 150, 'solar_zenith': 41, 'solar_azimuth': 67, + 'dni_extra': 1360, 'airmass': 1.324} + dni_adj, dhi_adj = ants2d._apply_sky_diffuse_model(**inputs, model=model) + # ensure that the adjusted values+isotropic yield the same poa_global as + # the target model + kwargs = inputs.copy() + kwargs.pop('dni') + kwargs.pop('dhi') + if model != 'perez': + kwargs.pop('airmass') + kwargs['surface_tilt'] = 20 + kwargs['surface_azimuth'] = 180 + adj = pvlib.irradiance.get_total_irradiance(dni=dni_adj, dhi=dhi_adj, + ghi=1000, # doesn't matter + model='isotropic', **kwargs) + func = {'perez': pvlib.irradiance.perez, + 'haydavies': pvlib.irradiance.haydavies}[model] + diffuse = func(dni=inputs['dni'], dhi=inputs['dhi'], **kwargs, + return_components=True) + aoi_proj = pvlib.irradiance.aoi_projection(kwargs['surface_tilt'], + kwargs['surface_azimuth'], + kwargs['solar_zenith'], + kwargs['solar_azimuth']) + poa_direct = inputs['dni'] * aoi_proj + diffuse['poa_circumsolar'] + poa_sky_diffuse = diffuse['poa_isotropic'] + poa_ground = 1000 * 0.25 * (1 - pvlib.tools.cosd(20)) / 2 + assert adj['poa_direct'] == pytest.approx(poa_direct, abs=1e-10) + assert adj['poa_sky_diffuse'] == pytest.approx(poa_sky_diffuse, abs=1e-10) + assert adj['poa_ground_diffuse'] == pytest.approx(poa_ground, abs=1e-10) + # sum of components, ignoring horizon brightening per ANTS-2D assumption + assert adj['poa_global'] == pytest.approx( + poa_direct + poa_sky_diffuse + poa_ground, abs=1e-10) + + +def test__apply_sky_diffuse_model_errors(): + with pytest.raises(ValueError, match='Must supply dni_extra'): + ants2d._apply_sky_diffuse_model(0, 0, 'haydavies', None, + None, None, None) + with pytest.raises(ValueError, match='Must supply dni_extra and airmass'): + ants2d._apply_sky_diffuse_model(0, 0, 'perez', None, + None, None, None) + with pytest.raises(ValueError, match='Invalid model: not_a_model'): + ants2d._apply_sky_diffuse_model(0, 0, 'not_a_model', None, + None, None, None) + + +def test__apply_ground_slope_cross_axis_slope(): + # use an outrageous cross_axis_slope, just for testing + inputs = { + 'height': 2, 'pitch': 4, 'gcr': 0.5, 'tracker_rotation': 50, + 'ghi': 600, 'dni': 900, 'dhi': 150, 'solar_zenith': 60, + 'solar_azimuth': 270, + } + outputs = ants2d._apply_ground_slope(axis_tilt=0, axis_azimuth=180, + cross_axis_slope=60, **inputs) + # cosd(60) gives a factor of 2 difference in various inputs + expected = {'height': inputs['height'] / 2, 'pitch': inputs['pitch'] * 2, + 'gcr': inputs['gcr'] / 2, 'tracker_rotation': -10, + # zenith aligns with cross_axis_slope: + 'ghi': inputs['dni'] + inputs['dhi']} + for (key, exp), actual in zip(expected.items(), outputs): + assert actual == pytest.approx(exp, abs=1e-10), key + + # now flip it around and check negative cross-axis slope + outputs = ants2d._apply_ground_slope(axis_tilt=0, axis_azimuth=180, + cross_axis_slope=-60, **inputs) + expected = {'height': inputs['height'] / 2, 'pitch': inputs['pitch'] * 2, + 'gcr': inputs['gcr'] / 2, 'tracker_rotation': 110, + 'ghi': inputs['dhi']} + for (key, exp), actual in zip(expected.items(), outputs): + assert actual == pytest.approx(exp, abs=1e-10), key + + +def test__apply_ground_slope_axis_tilt(): + # use an outrageous axis_tilt, just for testing + inputs = { + 'height': 2, 'pitch': 4, 'gcr': 0.5, 'tracker_rotation': 50, + 'ghi': 600, 'dni': 900, 'dhi': 150, 'solar_zenith': 60, + 'solar_azimuth': 180, + } + outputs = ants2d._apply_ground_slope(axis_tilt=60, axis_azimuth=180, + cross_axis_slope=0, **inputs) + expected = {'height': inputs['height'] / 2, 'pitch': inputs['pitch'], + 'gcr': inputs['gcr'], + 'tracker_rotation': inputs['tracker_rotation'], + # zenith aligns with axis_tilt: + 'ghi': inputs['dni'] + inputs['dhi']} + for (key, exp), actual in zip(expected.items(), outputs): + assert actual == pytest.approx(exp, abs=1e-10), key + + # now flip it around and check negative axis tilt + outputs = ants2d._apply_ground_slope(axis_tilt=-60, axis_azimuth=180, + cross_axis_slope=0, **inputs) + + expected = {'height': inputs['height'] / 2, 'pitch': inputs['pitch'], + 'gcr': inputs['gcr'], + 'tracker_rotation': inputs['tracker_rotation'], + 'ghi': inputs['dhi']} + for (key, exp), actual in zip(expected.items(), outputs): + assert actual == pytest.approx(exp, abs=1e-10), key + + +def test__apply_ground_slope_both(): + inputs = { + 'height': 2, 'pitch': 4, 'gcr': 0.5, 'tracker_rotation': 50, + 'ghi': 600, 'dni': 900, 'dhi': 150, 'solar_zenith': 15, + # the azimuth that results from the tilt/slope parameters below + 'solar_azimuth': 234.73561031724535, + } + outputs = ants2d._apply_ground_slope(axis_tilt=45, axis_azimuth=180, + cross_axis_slope=45, **inputs) + expected = {'height': inputs['height'] / 2, + 'pitch': inputs['pitch'] * 2**0.5, + 'gcr': inputs['gcr'] / 2**0.5, + 'tracker_rotation': inputs['tracker_rotation'] - 45, + # zenith aligns with axis_tilt: + 'ghi': inputs['dni'] / 2**0.5 + inputs['dhi']} + for (key, exp), actual in zip(expected.items(), outputs): + assert actual == pytest.approx(exp, abs=1e-10), key + + +def test__apply_ground_slope_zero(): + inputs = [2, 4, 0.5, 45, 600, 900, 150, 60, 210] + outputs = ants2d._apply_ground_slope(*inputs, axis_tilt=0, + axis_azimuth=180, + cross_axis_slope=0) + for input, output in zip(inputs, outputs): + assert pytest.approx(input, abs=1e-10) == output + + +@pytest.fixture +def ants_params(): + # parameters for get_irradiance + times = pd.date_range("2019-06-01 11:30", freq="h", periods=2) + inputs = { + 'tracker_rotation': [45, -45], + 'axis_azimuth': 180, + 'solar_zenith': [60, 60], + 'solar_azimuth': [225, 135], + 'gcr': 0.5, 'height': 2.5, 'pitch': 4.0, + 'ghi': [700, 700], + 'dni': [1000, 1000], + 'dhi': [200, 200], + 'albedo': 0.2, + 'dni_extra': [1360, 1360], + 'airmass': [2, 2], + } + for k, v in inputs.items(): + if isinstance(v, list): + inputs[k] = pd.Series(v, index=times) + return inputs + + +def test_get_irradiance_return_type(ants_params): + # verify pandas in -> pandas out, and shapes of numpy outputs + out = ants2d.get_irradiance(**ants_params, row_segments=1) + assert isinstance(out, pd.DataFrame) # DataFrame, since n_row_segments=1 + expected_keys = [ + 'poa_front', 'poa_front_direct', 'poa_front_diffuse', + 'poa_front_sky_diffuse', 'poa_front_ground_diffuse', + 'shaded_fraction_front', 'poa_back', 'poa_back_direct', + 'poa_back_diffuse', 'poa_back_sky_diffuse', 'poa_back_ground_diffuse', + 'shaded_fraction_back'] + assert set(out.columns) == set(expected_keys) + assert len(out) == 2 # 2 timestamps + + out = ants2d.get_irradiance(**ants_params, row_segments=3) + assert isinstance(out, dict) # dict, since row_segments>1 + assert set(out.keys()) == set(expected_keys) + for k, v in out.items(): + assert v.shape == (3, 2), k # 3 row segments, 2 timestamps + + +def test_get_irradiance_symmetry(ants_params): + # check symmetries for normal tracker + out = ants2d.get_irradiance(**ants_params, row_segments=1) + # symmetrical/mirrored inputs should produce equal outputs + pd.testing.assert_series_equal(out.iloc[0, :], out.iloc[1, :], + check_names=False) + + +@pytest.mark.parametrize('solar_zenith', [ + 60, # partial ground shading, no module shading + 80, # full ground shading, partial module shading +]) +@pytest.mark.parametrize('tracker_rotation', [+90, -90]) +def test_get_irradiance_vertical(ants_params, solar_zenith, tracker_rotation): + # check symmetries for vertical panels (tilt=90) + ants_params['solar_zenith'] = pd.Series(solar_zenith, + index=ants_params['ghi'].index) + ants_params['tracker_rotation'] = pd.Series(tracker_rotation, + index=ants_params['ghi'].index) + out = ants2d.get_irradiance(**ants_params, row_segments=1) + # inputs are symmetrical morning/afternoon, so morning front should equal + # afternoon back, and vice versa + front_keys = [ + 'poa_front', 'poa_front_direct', 'poa_front_diffuse', + 'poa_front_sky_diffuse', 'poa_front_ground_diffuse', + 'shaded_fraction_front'] + for front_key in front_keys: + back_key = front_key.replace("front", "back") + np.testing.assert_allclose(out.iloc[0][front_key], + out.iloc[1][back_key]) + np.testing.assert_allclose(out.iloc[1][front_key], + out.iloc[0][back_key]) + + # now with >1 row segment + out = ants2d.get_irradiance(**ants_params, row_segments=2) + lower_half = 0 + upper_half = 1 + morning = 0 + afternoon = 1 + for front_key in front_keys: + back_key = front_key.replace("front", "back") + np.testing.assert_allclose(out[front_key][lower_half, morning], + out[back_key][lower_half, afternoon]) + np.testing.assert_allclose(out[front_key][upper_half, morning], + out[back_key][upper_half, afternoon]) + np.testing.assert_allclose(out[back_key][lower_half, morning], + out[front_key][lower_half, afternoon]) + np.testing.assert_allclose(out[back_key][upper_half, morning], + out[front_key][upper_half, afternoon]) + + +def test_get_irradiance_limit(ants_params): + # check that diffuse components of front-side irradiance are lower + # than what get_total_irradiance predicts + surface_tilt = ants_params['tracker_rotation'].abs() + surface_azimuth = np.where(ants_params['tracker_rotation'] > 0, 270, 90) + irrad = pvlib.irradiance.get_total_irradiance( + surface_tilt, surface_azimuth, + ants_params['solar_zenith'], ants_params['solar_azimuth'], + ants_params['dni'], ants_params['ghi'], ants_params['dhi'], + albedo=ants_params['albedo'], model='isotropic') + + ants = ants2d.get_irradiance(**ants_params, row_segments=1, + model='isotropic') + # 15 W/m2 happens to be just below the difference (determined empirically) + diff_sky = irrad['poa_sky_diffuse'] - ants['poa_front_sky_diffuse'] + diff_grd = irrad['poa_ground_diffuse'] - ants['poa_front_ground_diffuse'] + assert all(diff_sky > 15) + assert all(diff_grd > 15) + + # but as pitch->infinity, front-side irradiance converges to + # output of get_total_irradiance + ants_params['pitch'] *= 1000 + ants_params['gcr'] /= 1000 + ants = ants2d.get_irradiance(**ants_params, row_segments=1, + model='isotropic') + + colmap = {'poa_front': 'poa_global', 'poa_front_direct': 'poa_direct', + 'poa_front_diffuse': 'poa_diffuse', + 'poa_front_sky_diffuse': 'poa_sky_diffuse', + 'poa_front_ground_diffuse': 'poa_ground_diffuse'} + ants_front = ants[list(colmap)].rename(columns=colmap) + pd.testing.assert_frame_equal(ants_front, irrad, atol=0.1) + + +@pytest.fixture +def ants_params_fixed(): + # parameters for get_irradiance, for a fixed-tilt system + inputs = { + 'tracker_rotation': 30, + 'axis_azimuth': 90, + 'solar_zenith': 60, + 'solar_azimuth': 175, + 'gcr': 0.6, 'height': 1.5, 'pitch': 3.5, + 'ghi': 700, + 'dni': 1000, + 'dhi': 200, + 'albedo': 0.2, + 'dni_extra': 1360, + 'airmass': 2, + } + return inputs + + +def test_get_irradiance_horizontal(ants_params_fixed): + # check that no issues pop up with tracker_rotation=solar_zenith=0 + ants_params_fixed['tracker_rotation'] = 0 + ants_params_fixed['solar_zenith'] = 0 + ants_params_fixed['solar_azimuth'] = 180 + zero = ants2d.get_irradiance(**ants_params_fixed) + + ants_params_fixed['tracker_rotation'] = 0.0001 + ants_params_fixed['solar_zenith'] = 0.0001 + ants_params_fixed['solar_azimuth'] = 180.0001 + pos_epsilon = ants2d.get_irradiance(**ants_params_fixed) + + ants_params_fixed['tracker_rotation'] = -0.0001 + ants_params_fixed['solar_zenith'] = -0.0001 + ants_params_fixed['solar_azimuth'] = 179.9999 + neg_epsilon = ants2d.get_irradiance(**ants_params_fixed) + + for key in zero: + assert_allclose(zero[key], pos_epsilon[key], atol=0.01) + assert_allclose(zero[key], neg_epsilon[key], atol=0.01) + + +def test_get_irradiance_direct_shading(ants_params_fixed): + # check that direct shading increases as sun approaches horizon + ants_params_fixed.pop('solar_zenith') + out60 = ants2d.get_irradiance(solar_zenith=60, **ants_params_fixed) + out80 = ants2d.get_irradiance(solar_zenith=80, **ants_params_fixed) + assert out80['poa_front_direct'] < out60['poa_front_direct'] + + +def test_get_irradiance_multiple_row_segments(ants_params_fixed): + # check that granular sims average to the same value as row_segments=1 + out4 = ants2d.get_irradiance(**ants_params_fixed, row_segments=4) + out2 = ants2d.get_irradiance(**ants_params_fixed, row_segments=2) + out1 = ants2d.get_irradiance(**ants_params_fixed, row_segments=1) + + for k in out4: + # check two bottom quarters average to the bottom half, and top + # two quarters average to the top half + np.testing.assert_allclose(np.mean(out4[k][0:2, 0]), out2[k][0, 0]) + np.testing.assert_allclose(np.mean(out4[k][2:4, 0]), out2[k][1, 0]) + + # check that two halves average to the whole + np.testing.assert_allclose(np.mean(out2[k][:, 0]), out1[k]) + + +def test_get_irradiance_custom_x0x1(ants_params_fixed): + # different ways of specifying the lower and upper halves + expected = ants2d.get_irradiance(**ants_params_fixed, row_segments=2) + actual = ants2d.get_irradiance(**ants_params_fixed, + row_segments=[(0.0, 0.5), (0.5, 1.0)]) + for key in expected: + np.testing.assert_allclose(expected[key], actual[key]) + + # specify only one part of the module + expected = ants2d.get_irradiance(**ants_params_fixed, row_segments=4) + actual = ants2d.get_irradiance(**ants_params_fixed, + row_segments=[(0.25, 0.5)]) + for key in expected: + np.testing.assert_allclose(expected[key][1], actual[key][0]) + + +def test_get_irradiance_slope(ants_params_fixed): + # check the slope affects direct & diffuse shading + flat = ants2d.get_irradiance(cross_axis_slope=0, **ants_params_fixed) + # negative slope with axis_azimuth=90 means sloping down to the north + tilt = ants2d.get_irradiance(cross_axis_slope=-10, **ants_params_fixed) + assert tilt['shaded_fraction_front'] > flat['shaded_fraction_front'] + assert tilt['poa_front_direct'] < flat['poa_front_direct'] + assert tilt['poa_front_sky_diffuse'] < flat['poa_front_sky_diffuse'] + + +def test_get_irradiance_nonuniform_albedo(): + # check that specifying albedo for each ground segment works + + # horizontal array, very close to the ground, with different albedo + # on left and right sides. check that ground-reflected irradiance at + # the edges of the module match the corresponding albedos + inputs = { + 'tracker_rotation': 0, + 'axis_azimuth': 180, + 'solar_zenith': 0, + 'solar_azimuth': 180, + 'gcr': 0.1, 'height': 0.05, 'pitch': 20, + 'ghi': 1000, + 'dni': 1000, + 'dhi': 0, + 'albedo': np.array([[0.5]*10 + [0.1]*10]).T, + 'model': 'isotropic' + } + out = ants2d.get_irradiance(ground_segments=20, + row_segments=10000, + max_rows=2, + **inputs) + # check far left and right segments, on the edge of the module. + # need a large row_segments so that these segments are very thin + left, right = out['poa_back_ground_diffuse'][[0, -1], 0] + # divide by two because ~half the visible ground is fully shaded + assert_allclose(left, 0.1 * 1000 / 2, rtol=0.002) + assert_allclose(right, 0.5 * 1000 / 2, rtol=0.002) + + +def test_get_irradiance_nonuniform_albedo_limit(): + # nonuniform albedo averages to uniform albedo, when sufficiently far away + base_inputs = { + 'tracker_rotation': 45, + 'axis_azimuth': 180, + 'solar_zenith': 10, + 'solar_azimuth': 215, + 'gcr': 0.5, 'height': 1000, 'pitch': 4, + 'ghi': 300, + 'dni': 0, # set dni to zero so that shadows don't confound results + 'dhi': 300, + 'ground_segments': 2, + 'max_rows': 10000, + 'model': 'isotropic', + } + out_uni = ants2d.get_irradiance(albedo=0.3, + **base_inputs) + out_non = ants2d.get_irradiance(albedo=np.array([[0.5, 0.1]]).T, + **base_inputs) + for key in out_non: + assert_allclose(out_non[key], out_uni[key], atol=1e-6) + + +@pytest.mark.parametrize('model,expected', [ + ('isotropic', { + 'poa_front': 1006.3548761345762, + 'poa_front_direct': 833.3333333333335, + 'poa_front_diffuse': 173.0215428012428, + 'poa_front_sky_diffuse': 172.27247024391784, + 'poa_front_ground_diffuse': 0.7490725573249604, + 'shaded_fraction_front': 0.035915234551783914, + 'poa_back': 23.626216052516494, + 'poa_back_direct': 0.0, + 'poa_back_diffuse': 23.626216052516494, + 'poa_back_sky_diffuse': 8.509173579096064, + 'poa_back_ground_diffuse': 15.11704247342043, + 'shaded_fraction_back': 0.035915234551784025}), + ('haydavies', { + 'poa_front': 1124.2311927022897, + 'poa_front_direct': 1078.4313725490197, + 'poa_front_diffuse': 45.79982015327015, + 'poa_front_sky_diffuse': 45.60153624103707, + 'poa_front_ground_diffuse': 0.19828391223307773, + 'shaded_fraction_front': 0.035915234551783914, + 'poa_back': 6.2539983668426, + 'poa_back_direct': 0.0, + 'poa_back_diffuse': 6.2539983668426, + 'poa_back_sky_diffuse': 2.252428300348958, + 'poa_back_ground_diffuse': 4.001570066493642, + 'shaded_fraction_back': 0.035915234551784025}), + ('perez', { + 'poa_front': 1060.3368384162613, + 'poa_front_direct': 945.5770264984124, + 'poa_front_diffuse': 114.75981191784896, + 'poa_front_sky_diffuse': 114.26297537137229, + 'poa_front_ground_diffuse': 0.4968365464766687, + 'shaded_fraction_front': 0.035915234551783914, + 'poa_back': 15.670534816764919, + 'poa_back_direct': 0.0, + 'poa_back_diffuse': 15.670534816764919, + 'poa_back_sky_diffuse': 5.643870374194696, + 'poa_back_ground_diffuse': 10.026664442570222, + 'shaded_fraction_back': 0.035915234551784025}) +]) +def test_get_irradiance_regression(model, expected, ants_params_fixed): + # values computed for typical but arbitrary inputs, to verify that output + # is stable over time + out = ants2d.get_irradiance(**ants_params_fixed, model=model) + for key in expected: + assert_allclose(out[key], expected[key], atol=1e-10) + + +def test_scalar_surface_angles(ants_params): + # scalar surface angles but Series irradiance inputs + ants_params['tracker_rotation'] = 30 + ants_params['axis_azimuth'] = 90 + out = ants2d.get_irradiance(**ants_params) + for key in out: + assert isinstance(out[key], pd.Series), key + assert len(out[key]) == 2, key + + # array inputs + ants_params = { + k: (v.values if isinstance(v, pd.Series) else v) + for k, v in ants_params.items() + } + out = ants2d.get_irradiance(**ants_params) + for key in out: + assert isinstance(out[key], np.ndarray), key + assert out[key].shape == (2,), key + + +def test_ground_components(): + # test ground incident irradiance values + inputs = { + 'solar_zenith': 60, + 'solar_azimuth': 90, + 'gcr': 0.5, 'height': 2.5, 'pitch': 4, + 'ghi': 800, + 'dni': 1000, + 'dhi': 300, + 'albedo': 0.2, + 'ground_segments': 1, + 'model': 'isotropic', + } + + # sun is edge-on to modules: ground is fully illuminated + _, out = ants2d.get_irradiance(tracker_rotation=30, axis_azimuth=180, + **inputs, + return_ground_components=True) + assert_allclose(out['ground_direct'], 500) # dni * cos(zenith) + assert out['ground_diffuse'] < 300 # some dhi blocked by rows + + # sun is normal to modules: ground is fully shaded + _, out = ants2d.get_irradiance(tracker_rotation=-60, axis_azimuth=180, + **inputs, + return_ground_components=True) + assert_allclose(out['ground_direct'], 0, atol=1e-10) + assert out['ground_diffuse'] < 300 + + # flat array: ground_direct -> dni * cos(zenith) * gcr + _, out = ants2d.get_irradiance(tracker_rotation=0, axis_azimuth=180, + **inputs, + return_ground_components=True, + max_rows=500) + assert_allclose(out['ground_direct'], 250) # gcr of 0.5 + assert_allclose(out['ground_diffuse'], 150, atol=1e-3) + + # flat array, four segments: perfect direct shading, diffuse symmetry + inputs['ground_segments'] = 4 + inputs['solar_zenith'] = 0 + _, out = ants2d.get_irradiance(tracker_rotation=0, axis_azimuth=180, + **inputs, + return_ground_components=True) + assert_allclose(out['ground_direct'], [0, 500, 500, 0]) + diffuse = out['ground_diffuse'] + assert_allclose(diffuse[0], diffuse[3], atol=0.1) + assert_allclose(diffuse[1], diffuse[2], atol=0.1) + + # flat array, many rows, very high: ground_diffuse -> dhi * gcr + inputs['height'] = 200 + inputs['ground_segments'] = 4 + _, out = ants2d.get_irradiance(tracker_rotation=0, axis_azimuth=180, + **inputs, + max_rows=5000, + return_ground_components=True) + assert_allclose(out['ground_diffuse'], 150, atol=0.01) + + +def test_ground_components_types(ants_params, ants_params_fixed): + # test second return value type/shape when return_ground_components=True + + # scalar inputs, single ground segment + _, out = ants2d.get_irradiance(**ants_params_fixed, ground_segments=1, + return_ground_components=True) + assert isinstance(out, dict) + assert set(out.keys()) == {'ground_direct', 'ground_diffuse'} + for key, value in out.items(): + assert isinstance(value, float), key + + # scalar inputs, multiple ground segments + _, out = ants2d.get_irradiance(**ants_params_fixed, ground_segments=10, + return_ground_components=True) + assert isinstance(out, dict) + assert set(out.keys()) == {'ground_direct', 'ground_diffuse'} + for key, value in out.items(): + assert isinstance(value, np.ndarray), key + assert value.shape == (10,), key + + # series inputs, single ground segment + _, out = ants2d.get_irradiance(**ants_params, ground_segments=1, + return_ground_components=True) + assert isinstance(out, pd.DataFrame) + assert set(out.columns) == {'ground_direct', 'ground_diffuse'} + assert len(out) == 2 + + # series inputs, multiple ground segments + _, out = ants2d.get_irradiance(**ants_params, ground_segments=10, + return_ground_components=True) + assert isinstance(out, dict) + assert set(out.keys()) == {'ground_direct', 'ground_diffuse'} + for key, value in out.items(): + assert isinstance(value, np.ndarray), key + assert value.shape == (10, 2), key diff --git a/tests/bifacial/test_infinite_sheds.py b/tests/bifacial/test_infinite_sheds.py index 5459b9d575..d91ef322e0 100644 --- a/tests/bifacial/test_infinite_sheds.py +++ b/tests/bifacial/test_infinite_sheds.py @@ -9,6 +9,8 @@ import pytest +from pvlib._deprecation import pvlibDeprecationWarning + @pytest.fixture def test_system(): @@ -94,10 +96,12 @@ def test_get_irradiance_poa(): albedo = 0 iam = 1.0 npoints = 100 - res = infinite_sheds.get_irradiance_poa( - surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, - gcr, height, pitch, ghi, dhi, dni, - albedo, iam=iam, npoints=npoints) + match = "`npoints` and `vectorize` parameters have no effect" + with pytest.warns(pvlibDeprecationWarning, match=match): + res = infinite_sheds.get_irradiance_poa( + surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, + gcr, height, pitch, ghi, dhi, dni, + albedo, iam=iam, npoints=npoints) expected_diffuse = np.array([300.]) expected_direct = np.array([700.]) expected_global = expected_diffuse + expected_direct @@ -120,10 +124,12 @@ def test_get_irradiance_poa(): expected_global = expected_diffuse + expected_direct expected_shaded_fraction = np.array( [0., 0., 0., 0.]) - res = infinite_sheds.get_irradiance_poa( - surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, - gcr, height, pitch, ghi, dhi, dni, - albedo, iam=iam, npoints=npoints) + match = "`npoints` and `vectorize` parameters have no effect" + with pytest.warns(pvlibDeprecationWarning, match=match): + res = infinite_sheds.get_irradiance_poa( + surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, + gcr, height, pitch, ghi, dhi, dni, + albedo, iam=iam, npoints=npoints) assert np.allclose(res['poa_global'], expected_global) assert np.allclose(res['poa_diffuse'], expected_diffuse) assert np.allclose(res['poa_direct'], expected_direct) @@ -143,10 +149,12 @@ def test_get_irradiance_poa(): expected_shaded_fraction = pd.Series( data=expected_shaded_fraction, index=surface_tilt.index) expected_shaded_fraction.name = 'shaded_fraction' # to match output Series - res = infinite_sheds.get_irradiance_poa( - surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, - gcr, height, pitch, ghi, dhi, dni, - albedo, iam=iam, npoints=npoints) + match = "`npoints` and `vectorize` parameters have no effect" + with pytest.warns(pvlibDeprecationWarning, match=match): + res = infinite_sheds.get_irradiance_poa( + surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, + gcr, height, pitch, ghi, dhi, dni, + albedo, iam=iam, npoints=npoints) assert isinstance(res, pd.DataFrame) assert_series_equal(res['poa_global'], expected_global) assert_series_equal(res['shaded_fraction'], expected_shaded_fraction) @@ -180,11 +188,13 @@ def test_get_irradiance(vectorize): iam_front = 1.0 iam_back = 1.0 npoints = 100 - result = infinite_sheds.get_irradiance( - surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, - gcr, height, pitch, ghi, dhi, dni, albedo, iam_front, iam_back, - bifaciality=0.8, shade_factor=-0.02, transmission_factor=0, - npoints=npoints, vectorize=vectorize) + match = "`npoints` and `vectorize` parameters have no effect" + with pytest.warns(pvlibDeprecationWarning, match=match): + result = infinite_sheds.get_irradiance( + surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, + gcr, height, pitch, ghi, dhi, dni, albedo, iam_front, iam_back, + bifaciality=0.8, shade_factor=-0.02, transmission_factor=0, + npoints=npoints, vectorize=vectorize) expected_front_diffuse = np.array([300.]) expected_front_direct = np.array([700.]) expected_front_global = expected_front_diffuse + expected_front_direct @@ -204,15 +214,17 @@ def test_get_irradiance(vectorize): dni = pd.Series([700., 0., 0., 700.], index=ghi.index) solar_zenith = pd.Series([0., 0., 0., 135.], index=ghi.index) surface_tilt = pd.Series([0., 0., 90., 0.], index=ghi.index) - result = infinite_sheds.get_irradiance( - surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, - gcr, height, pitch, ghi, dhi, dni, albedo, iam_front, iam_back, - bifaciality=0.8, shade_factor=-0.02, transmission_factor=0, - npoints=npoints, vectorize=vectorize) - result_front = infinite_sheds.get_irradiance_poa( - surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, - gcr, height, pitch, ghi, dhi, dni, - albedo, iam=iam_front, vectorize=vectorize) + match = "`npoints` and `vectorize` parameters have no effect" + with pytest.warns(pvlibDeprecationWarning, match=match): + result = infinite_sheds.get_irradiance( + surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, + gcr, height, pitch, ghi, dhi, dni, albedo, iam_front, iam_back, + bifaciality=0.8, shade_factor=-0.02, transmission_factor=0, + npoints=npoints, vectorize=vectorize) + result_front = infinite_sheds.get_irradiance_poa( + surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, + gcr, height, pitch, ghi, dhi, dni, + albedo, iam=iam_front, vectorize=vectorize) assert isinstance(result, pd.DataFrame) expected_poa_global = pd.Series( [1000., 500., result_front['poa_global'][2] * (1 + 0.8 * 0.98), @@ -234,7 +246,7 @@ def test_get_irradiance_limiting_gcr(): surface_azimuth = 180. gcr = 0.00001 height = 1. - pitch = 100. + pitch = 1 / gcr ghi = 1000. dhi = 300. dni = 700. @@ -242,11 +254,13 @@ def test_get_irradiance_limiting_gcr(): iam_front = 1.0 iam_back = 1.0 npoints = 100 - result = infinite_sheds.get_irradiance( - surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, - gcr, height, pitch, ghi, dhi, dni, albedo, iam_front, iam_back, - bifaciality=1., shade_factor=-0.00, transmission_factor=0., - npoints=npoints) + match = "`npoints` and `vectorize` parameters have no effect" + with pytest.warns(pvlibDeprecationWarning, match=match): + result = infinite_sheds.get_irradiance( + surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, + gcr, height, pitch, ghi, dhi, dni, albedo, iam_front, iam_back, + bifaciality=1., shade_factor=-0.00, transmission_factor=0., + npoints=npoints) expected_ground_diffuse = np.array([500.]) expected_sky_diffuse = np.array([150.]) expected_direct = np.array([0.]) @@ -292,11 +306,13 @@ def test_get_irradiance_with_haydavies(): iam_front = 1.0 iam_back = 1.0 npoints = 100 - result = infinite_sheds.get_irradiance( - surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, - gcr, height, pitch, ghi, dhi, dni, albedo, model, dni_extra, - iam_front, iam_back, bifaciality=0.8, shade_factor=-0.02, - transmission_factor=0, npoints=npoints) + match = "`npoints` and `vectorize` parameters have no effect" + with pytest.warns(pvlibDeprecationWarning, match=match): + result = infinite_sheds.get_irradiance( + surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, + gcr, height, pitch, ghi, dhi, dni, albedo, model, dni_extra, + iam_front, iam_back, bifaciality=0.8, shade_factor=-0.02, + transmission_factor=0, npoints=npoints) expected_front_diffuse = np.array([151.38]) expected_front_direct = np.array([848.62]) expected_front_global = expected_front_diffuse + expected_front_direct diff --git a/tests/bifacial/test_utils.py b/tests/bifacial/test_utils.py index e24af42b3e..c0a6830167 100644 --- a/tests/bifacial/test_utils.py +++ b/tests/bifacial/test_utils.py @@ -8,6 +8,8 @@ from pvlib.tools import cosd from scipy.integrate import trapezoid +from pvlib._deprecation import pvlibDeprecationWarning + @pytest.fixture def test_system_fixed_tilt(): @@ -53,33 +55,42 @@ def test__solar_projection_tangent(): @pytest.mark.parametrize( - "gcr,surface_tilt,surface_azimuth,solar_zenith,solar_azimuth,expected", - [(0.5, 0., 180., 0., 180., 0.5), - (1.0, 0., 180., 0., 180., 0.0), - (1.0, 90., 180., 0., 180., 1.0), - (0.5, 45., 180., 45., 270., 1.0 - np.sqrt(2) / 4), - (0.5, 45., 180., 90., 180., 0.), - (np.sqrt(2) / 2, 45, 180, 0, 180, 0.5), - (np.sqrt(2) / 2, 45, 180, 45, 180, 0.0), - (np.sqrt(2) / 2, 45, 180, 45, 90, 0.5), - (np.sqrt(2) / 2, 45, 180, 45, 0, 1.0), - (np.sqrt(2) / 2, 45, 180, 45, 135, 0.5 * (1 - np.sqrt(2) / 2)), + "gcr,surface_tilt,phi,expected", + [(0.5, 0., 0., 0.5), + (1.0, 0., 0., 0.0), + (1.0, 90., 0., 1.0), + (0.5, 45., 0., 1.0 - np.sqrt(2) / 4), + (0.5, 45., 90., 0.), + (np.sqrt(2) / 2, 45, 0, 0.5), + (np.sqrt(2) / 2, 45, 45, 0.0), + (np.sqrt(2) / 2, 45, 0, 0.5), + (np.sqrt(2) / 2, 45, -45, 1.0), + (np.sqrt(2) / 2, 45, 35.264389682754654, 0.5 * (1 - np.sqrt(2) / 2)), ]) -def test__unshaded_ground_fraction( - surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, gcr, - expected): +def test__unshaded_ground_fraction(surface_tilt, phi, gcr, expected): # frontside, same for both sides - f_sky_beam_f = utils._unshaded_ground_fraction( - surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, gcr) + f_sky_beam_f = utils._unshaded_ground_fraction(surface_tilt, phi, gcr) assert np.allclose(f_sky_beam_f, expected) # backside, should be the same as frontside - f_sky_beam_b = utils._unshaded_ground_fraction( - 180. - surface_tilt, surface_azimuth - 180., solar_zenith, - solar_azimuth, gcr) + f_sky_beam_b = utils._unshaded_ground_fraction(surface_tilt - 180., phi, + gcr) assert np.allclose(f_sky_beam_b, expected) -def test__vf_ground_sky_2d(test_system_fixed_tilt): +def test__unshaded_ground_fraction_horizontal(): + # test that zero angles don't mess things up. specifically, check + # for continuity with small positive and negative angles + params = dict(gcr=0.5, pitch=4, height=2.5, + g0=[0, 0.25, 0.5, 0.75], g1=[0.25, 0.5, 0.75, 1]) + zero = utils._unshaded_ground_fraction(0.0, 0.0, **params) + pos = utils._unshaded_ground_fraction(0.01, 0.01, **params) + neg = utils._unshaded_ground_fraction(-0.01, -0.01, **params) + np.testing.assert_allclose(pos, neg, atol=0.01) + np.testing.assert_allclose(pos, zero, atol=0.01) + np.testing.assert_allclose(neg, zero, atol=0.01) + + +def test_vf_ground_sky_2d(test_system_fixed_tilt): # vector input ts, pts, vfs_gnd_sky = test_system_fixed_tilt vfs = utils.vf_ground_sky_2d(ts['rotation'], ts['gcr'], pts, @@ -97,9 +108,11 @@ def test_vf_ground_sky_2d_integ(test_system_fixed_tilt, vectorize): # pass rotation here since max_rows=1 for the hand-solved case in # the fixture test_system, which means the ground-to-sky view factor # isn't summed over enough rows for symmetry to hold. - vf_integ = utils.vf_ground_sky_2d_integ( - ts['rotation'], ts['gcr'], ts['height'], ts['pitch'], - max_rows=1, npoints=3, vectorize=vectorize) + match = '`npoints` and `vectorize` parameters have no effect' + with pytest.warns(pvlibDeprecationWarning, match=match): + vf_integ = utils.vf_ground_sky_2d_integ( + ts['rotation'], ts['gcr'], ts['height'], ts['pitch'], + max_rows=1, npoints=3, vectorize=vectorize) expected_vf_integ = trapezoid(vfs_gnd_sky, pts, axis=0) assert np.isclose(vf_integ, expected_vf_integ, rtol=0.1) @@ -147,6 +160,19 @@ def test_vf_row_sky_2d_integ(test_system_fixed_tilt): assert np.allclose(vf, y1, rtol=1e-3) +def test_vf_ground_sky_2d_integ_horizontal_large_max_rows(): + # with large max_rows, rows far out towards the horizon are considered. + # obstructed string lengths are then calculated using angles very close + # to zero. however, numerical roundoff requires some amount of slop being + # allowed in the comparison. this case previously failed when the slop + # allowance was too large and prevented accurate resolution of those + # far-away rows. + inputs = dict(tracker_rotation=0, gcr=0.518, height=1.5, pitch=3.5, + g0=0, g1=0.05, max_rows=1000) + vf_gnd_sky = utils.vf_ground_sky_2d_integ(**inputs) + assert np.max(vf_gnd_sky) < 1 + + def test_vf_row_ground_2d(test_system_fixed_tilt): ts, _, _ = test_system_fixed_tilt # with float input, fx at bottom of row @@ -161,20 +187,22 @@ def test_vf_row_ground_2d(test_system_fixed_tilt): assert np.allclose(vf, expected) -def test_vf_ground_2d_integ(test_system_fixed_tilt): +def test_vf_row_ground_2d_integ(test_system_fixed_tilt): ts, _, _ = test_system_fixed_tilt # with float input, check end position with np.errstate(invalid='ignore'): vf = utils.vf_row_ground_2d_integ(ts['surface_tilt'], ts['gcr'], - 0., 0.) + ts['height'], ts['pitch'], + x0=0., x1=0.00001, max_rows=1000) expected = utils.vf_row_ground_2d(ts['surface_tilt'], ts['gcr'], 0.) - assert np.isclose(vf, expected) + assert np.isclose(vf, expected, atol=1e-4) # with array input fx0 = np.array([0., 0.5]) - fx1 = np.array([0., 0.8]) + fx1 = np.array([0.00001, 0.8]) with np.errstate(invalid='ignore'): vf = utils.vf_row_ground_2d_integ(ts['surface_tilt'], ts['gcr'], - fx0, fx1) + ts['height'], ts['pitch'], + x0=fx0, x1=fx1, max_rows=1000) phi = ground_angle(ts['surface_tilt'], ts['gcr'], fx0[0]) y0 = 0.5 * (1 - cosd(phi - ts['surface_tilt'])) x = np.arange(fx0[1], fx1[1], 1e-4) @@ -184,9 +212,26 @@ def test_vf_ground_2d_integ(test_system_fixed_tilt): expected = np.array([y0, y1]) assert np.allclose(vf, expected, rtol=1e-2) # with defaults (0, 1) - vf = utils.vf_row_ground_2d_integ(ts['surface_tilt'], ts['gcr'], 0, 1) + vf = utils.vf_row_ground_2d_integ(ts['surface_tilt'], ts['gcr'], + ts['height'], ts['pitch'], x0=0, x1=1, + max_rows=1000) x = np.arange(0, 1, 1e-4) phi_y = ground_angle(ts['surface_tilt'], ts['gcr'], x) y = 0.5 * (1 - cosd(phi_y - ts['surface_tilt'])) y1 = trapezoid(y, x) / (1 - 0) assert np.allclose(vf, y1, rtol=1e-2) + + +def test_vf_row_ground_2d_integ_upsidedown(): + # horizontal, rear-side + inputs = {'surface_tilt': np.array([-180]), + 'gcr': 0.6, 'height': 1.5, 'pitch': 3.5, + 'x0': np.array([0.]), 'x1': np.array([1.]), + 'g0': np.array([0., 0.5]), + 'g1': np.array([0.5, 1.]), + 'max_rows': 7} + vf_row_ground = utils.vf_row_ground_2d_integ(**inputs) + assert all(vf_row_ground > 1e-3) + inputs['surface_tilt'] = +180 + vf_row_ground = utils.vf_row_ground_2d_integ(**inputs) + assert all(vf_row_ground > 1e-3)