pvfactors limited implementation for bifacial calculations

pvlib/bifacial.py

@@ -0,0 +1,135 @@
+"""
+The ``bifacial`` module contains functions for modeling back surface
+plane-of-array irradiance under various conditions.
+"""
+
+import pandas as pd
+from pvlib.tools import enforce_numpy_arrays
+
+
+@enforce_numpy_arrays
+def pvfactors_timeseries(
+        solar_azimuth, solar_zenith, surface_azimuth, surface_tilt,
+        timestamps, dni, dhi, gcr, pvrow_height, pvrow_width, albedo,
+        n_pvrows=3, index_observed_pvrow=1,
+        rho_front_pvrow=0.03, rho_back_pvrow=0.05,
+        horizon_band_angle=15.,
+        run_parallel_calculations=True, n_workers_for_parallel_calcs=None):
+    """
+    Calculate front and back surface plane-of-array irradiance on
+    a fixed tilt or single-axis tracker PV array configuration, and using
+    the open-source "pvfactors" package.
+    Please refer to pvfactors online documentation for more details:
+    https://sunpower.github.io/pvfactors/
+
+    Inputs
+    ------
+    solar_azimuth: numeric
+        Sun's azimuth angles using pvlib's azimuth convention (deg)
+    solar_zenith: numeric
+        Sun's zenith angles (deg)
+    surface_azimuth: numeric
+        Azimuth angle of the front surface of the PV modules, using pvlib's
+        convention (deg)
+    surface_tilt: numeric
+        Tilt angle of the PV modules, going from 0 to 180 (deg)
+    timestamps: array of :class:datetime.datetime objects
+        List of simulation timestamps
+    dni: numeric
+        Direct normal irradiance (W/m2)
+    dhi: numeric
+        Diffuse horizontal irradiance (W/m2)
+    gcr: float
+        Ground coverage ratio of the pv array
+    pvrow_height: float
+        Height of the pv rows, measured at their center (m)
+    pvrow_width: float
+        Width of the pv rows in the considered 2D plane (m)
+    albedo: float
+        Ground albedo
+    n_pvrows: int, default 3
+        Number of PV rows to consider in the PV array
+    index_observed_pvrow: int, default 1
+        Index of the PV row whose incident irradiance will be returned. Indices
+        of PV rows go from 0 to n_pvrows-1.
+    rho_front_pvrow: float, default 0.03
+        Front surface reflectivity of PV rows
+    rho_back_pvrow: float, default 0.05
+        Back surface reflectivity of PV rows
+    horizon_band_angle: float, default 15
+        Elevation angle of the sky dome's diffuse horizon band (deg)
+    run_parallel_calculations: bool, default True
+        pvfactors is capable of using multiprocessing. Use this flag to decide
+        to run calculations in parallel (recommended) or not.
+    n_workers_for_parallel_calcs: int, default None
+        Number of workers to use in the case of parallel calculations. The
+        default value of 'None' will lead to using a value equal to the number
+        of CPU's on the machine running the model.
+
+    Returns
+    -------
+    front_poa_irradiance: numeric
+        Calculated incident irradiance on the front surface of the PV modules
+        (W/m2)
+    back_poa_irradiance: numeric
+        Calculated incident irradiance on the back surface of the PV modules
+        (W/m2)
+    df_registries: pandas DataFrame
+        DataFrame containing all the detailed outputs and elements calculated
+        for every timestamp of the simulation. Please refer to pvfactors
+        documentation for more details
+
+    References
+    ----------
+    .. [1] Anoma, Marc Abou, et al. "View Factor Model and Validation for
+        Bifacial PV and Diffuse Shade on Single-Axis Trackers." 44th IEEE
+        Photovoltaic Specialist Conference. 2017.
+    """
+
+    # Import pvfactors functions for timeseries calculations.
+    from pvfactors.timeseries import (calculate_radiosities_parallel_perez,
+                                      calculate_radiosities_serially_perez,
+                                      get_average_pvrow_outputs)
+    idx_slice = pd.IndexSlice
+
+    # Build up pv array configuration parameters
+    pvarray_parameters = {
+        'n_pvrows': n_pvrows,
+        'pvrow_height': pvrow_height,
+        'pvrow_width': pvrow_width,
+        'surface_azimuth': surface_azimuth[0],  # not necessary
+        'surface_tilt': surface_tilt[0],        # not necessary
+        'gcr': gcr,
+        'solar_zenith': solar_zenith[0],        # not necessary
+        'solar_azimuth': solar_azimuth[0],      # not necessary
+        'rho_ground': albedo,
+        'rho_front_pvrow': rho_front_pvrow,
+        'rho_back_pvrow': rho_back_pvrow,
+        'horizon_band_angle': horizon_band_angle
+    }
+
+    # Run pvfactors calculations: either in parallel or serially
+    if run_parallel_calculations:
+        df_registries, df_custom_perez = calculate_radiosities_parallel_perez(
+            pvarray_parameters, timestamps, solar_zenith, solar_azimuth,
+            surface_tilt, surface_azimuth, dni, dhi,
+            n_processes=n_workers_for_parallel_calcs)
+    else:
+        inputs = (pvarray_parameters, timestamps, solar_zenith, solar_azimuth,
+                  surface_tilt, surface_azimuth, dni, dhi)
+        df_registries, df_custom_perez = calculate_radiosities_serially_perez(
+            inputs)
+
+    # Get the average surface outputs
+    df_outputs = get_average_pvrow_outputs(df_registries,
+                                           values=['qinc'],
+                                           include_shading=True)
+
+    # Select the calculated outputs from the pvrow to observe
+    ipoa_front = df_outputs.loc[:, idx_slice[index_observed_pvrow,
+                                             'front', 'qinc']].values
+
+    ipoa_back = df_outputs.loc[:, idx_slice[index_observed_pvrow,
+                                            'back', 'qinc']].values
+
+    return ipoa_front, ipoa_back, df_registries

pvlib/test/conftest.py

@@ -69,13 +69,15 @@ def inner():
 def pandas_0_17():
     return parse_version(pd.__version__) >= parse_version('0.17.0')
 
+
 needs_pandas_0_17 = pytest.mark.skipif(
     not pandas_0_17(), reason='requires pandas 0.17 or greater')
 
 
 def numpy_1_10():
     return parse_version(np.__version__) >= parse_version('1.10.0')
 
+
 needs_numpy_1_10 = pytest.mark.skipif(
     not numpy_1_10(), reason='requires numpy 1.10 or greater')
 
@@ -92,8 +94,10 @@ def has_spa_c():
     else:
         return True
 
+
 requires_spa_c = pytest.mark.skipif(not has_spa_c(), reason="requires spa_c")
 
+
 def has_numba():
     try:
         import numba
@@ -106,6 +110,7 @@ def has_numba():
         else:
             return True
 
+
 requires_numba = pytest.mark.skipif(not has_numba(), reason="requires numba")
 
 try:
@@ -125,3 +130,30 @@ def has_numba():
 
 requires_netCDF4 = pytest.mark.skipif(not has_netCDF4,
                                       reason='requires netCDF4')
+
+try:
+    import pvfactors
+    has_pvfactors = True
+except ImportError:
+    has_pvfactors = False
+
+requires_pvfactors = pytest.mark.skipif(not has_pvfactors,
+                                        reason='requires pvfactors')
+
+try:
+    import future
+    has_future = True
+except ImportError:
+    has_future = False
+
+requires_future = pytest.mark.skipif(not has_future,
+                                     reason='requires future')
+
+try:
+    import shapely
+    has_shapely = True
+except ImportError:
+    has_shapely = False
+
+requires_shapely = pytest.mark.skipif(not has_shapely,
+                                      reason='requires shapely')

pvlib/test/test_bifacial.py

@@ -0,0 +1,66 @@
+import numpy as np
+from datetime import datetime
+from pvlib.bifacial import pvfactors_timeseries
+from conftest import (requires_pvfactors, requires_future, requires_shapely,
+                      requires_scipy)
+
+
+@requires_scipy
+@requires_shapely
+@requires_future
+@requires_pvfactors
+def test_pvfactors_timeseries():
+    """ Test that pvfactors is functional, using the TLDR section inputs of the
+    package github repo README.md file"""
+
+    # Create some inputs
+    timestamps = np.array([datetime(2017, 8, 31, 11),
+                           datetime(2017, 8, 31, 12)])
+    solar_zenith = np.array([20., 10.])
+    solar_azimuth = np.array([110., 140.])
+    surface_tilt = np.array([10., 0.])
+    surface_azimuth = np.array([90., 90.])
+    dni = np.array([1000., 300.])
+    dhi = np.array([50., 500.])
+    gcr = 0.4
+    pvrow_height = 1.75
+    pvrow_width = 2.44
+    albedo = 0.2
+    n_pvrows = 3
+    index_observed_pvrow = 1
+    rho_front_pvrow = 0.03
+    rho_back_pvrow = 0.05
+    horizon_band_angle = 15.
+
+    # Expected values
+    expected_ipoa_front = [1034.96216923, 795.4423259]
+    expected_ipoa_back = [92.11871485, 70.39404124]
+    tolerance = 1e-6
+
+    # Test serial calculations
+    ipoa_front, ipoa_back, _ = pvfactors_timeseries(
+        solar_azimuth, solar_zenith, surface_azimuth, surface_tilt,
+        timestamps, dni, dhi, gcr, pvrow_height, pvrow_width, albedo,
+        n_pvrows=n_pvrows, index_observed_pvrow=index_observed_pvrow,
+        rho_front_pvrow=rho_front_pvrow, rho_back_pvrow=rho_back_pvrow,
+        horizon_band_angle=horizon_band_angle,
+        run_parallel_calculations=False, n_workers_for_parallel_calcs=None)
+
+    np.testing.assert_allclose(ipoa_front, expected_ipoa_front,
+                               atol=0, rtol=tolerance)
+    np.testing.assert_allclose(ipoa_back, expected_ipoa_back,
+                               atol=0, rtol=tolerance)
+
+    # Run calculations in parallel
+    ipoa_front, ipoa_back, _ = pvfactors_timeseries(
+        solar_azimuth, solar_zenith, surface_azimuth, surface_tilt,
+        timestamps, dni, dhi, gcr, pvrow_height, pvrow_width, albedo,
+        n_pvrows=n_pvrows, index_observed_pvrow=index_observed_pvrow,
+        rho_front_pvrow=rho_front_pvrow, rho_back_pvrow=rho_back_pvrow,
+        horizon_band_angle=horizon_band_angle,
+        run_parallel_calculations=True, n_workers_for_parallel_calcs=None)
+
+    np.testing.assert_allclose(ipoa_front, expected_ipoa_front,
+                               atol=0, rtol=tolerance)
+    np.testing.assert_allclose(ipoa_back, expected_ipoa_back,
+                               atol=0, rtol=tolerance)

pvlib/test/test_tools.py

@@ -1,6 +1,8 @@
 import pytest
 
 from pvlib import tools
+import numpy as np
+import pandas as pd
 
 @pytest.mark.parametrize('keys, input_dict, expected', [
     (['a', 'b'], {'a': 1, 'b': 2, 'c': 3}, {'a': 1, 'b': 2}),
@@ -11,3 +13,45 @@
 def test_build_kwargs(keys, input_dict, expected):
     kwargs = tools._build_kwargs(keys, input_dict)
     assert kwargs == expected
+
+
+def test_enforce_numpy_arrays():
+    """ Check that when pandas objects are included in the inputs, the
+    wrapped function receives numpy arrays and the outputs are turned into
+    pandas series """
+
+    # Create a test function and decorate it
+    @tools.enforce_numpy_arrays
+    def fun_function(a, b, c, d, kwarg_1=None, kwarg_2=None):
+        assert isinstance(a, np.ndarray)
+        assert isinstance(b, np.ndarray)
+        assert isinstance(c, np.ndarray)
+        assert isinstance(d, str)
+        assert isinstance(kwarg_1, np.ndarray)
+        assert isinstance(kwarg_2, np.ndarray)
+        out_1 = np.array([1, 2])
+        out_2 = pd.Series([3, 4])
+        out_3 = 5.
+        return out_1, out_2, out_3
+
+    # Check with no pandas inputs
+    a = b = c = np.array([1, 2])
+    d = 'string'
+    kwarg_1 = np.array([1, 2])
+    kwarg_2 = np.array([1, 2])
+    out_1, out_2, out_3 = fun_function(a, b, c, d,
+                                       kwarg_1=kwarg_1, kwarg_2=kwarg_2)
+    assert isinstance(out_1, np.ndarray)
+    assert isinstance(out_2, pd.Series)
+    assert isinstance(out_3, float)
+
+    # Check with some pandas inputs in both args and kwargs
+    b = pd.Series([1, 2])
+    c = pd.DataFrame([1, 2], columns=['e'], index=range(2))
+    kwarg_1 = pd.Series([1, 2])
+    kwarg_2 = pd.DataFrame([1, 2], columns=['kwarg_2'], index=range(2))
+    out_1, out_2, out_3 = fun_function(a, b, c, d,
+                                       kwarg_1=kwarg_1, kwarg_2=kwarg_2)
+    assert isinstance(out_1, pd.Series)
+    assert isinstance(out_2, pd.Series)
+    assert isinstance(out_3, float)

pvlib/tools.py

@@ -8,6 +8,7 @@
 import numpy as np
 import pandas as pd
 import pytz
+from functools import wraps
 
 
 def cosd(angle):
@@ -425,3 +426,58 @@ def _golden_sect_DataFrame(params, VL, VH, func):
             raise Exception("EXCEPTION:iterations exceeded maximum (50)")
 
     return func(df, 'V1'), df['V1']
+
+
+def enforce_numpy_arrays(f):
+    """ This decorator function will convert any inputted pandas dataframe or
+    series to a numpy array, and make sure to convert outtputed numpy arrays
+    to pandas series if any input was a pandas structure """
+
+    @wraps(f)
+    def wrapper(*args, **kwargs):
+        # Convert all inputs to numpy arrays if pandas series or dataframe
+        is_return_type_pandas = False
+
+        new_args = []
+        for arg in args:
+            if isinstance(arg, pd.Series):
+                is_return_type_pandas = True
+                new_arg = arg.values
+                new_args.append(new_arg)
+            elif isinstance(arg, pd.DataFrame):
+                is_return_type_pandas = True
+                new_arg = arg.values.ravel()  # make sure that 1D
+                new_args.append(new_arg)
+            else:
+                new_args.append(arg)
+        new_args = tuple(new_args)
+
+        new_kwargs = {}
+        for key, val in kwargs.items():
+            if isinstance(val, pd.DataFrame) or isinstance(val, pd.Series):
+                is_return_type_pandas = True
+                new_val = val.values
+                new_kwargs[key] = new_val
+            else:
+                new_kwargs[key] = val
+
+        # Run function
+        outputs = f(*new_args, **new_kwargs)
+
+        # If inputs were provided as pandas series or dataframe, convert
+        # outputted numpy arrays to pandas series
+        new_outputs = []
+        if is_return_type_pandas:
+            for output in outputs:
+                if isinstance(output, np.ndarray):
+                    new_output = pd.Series(output)
+                    new_outputs.append(new_output)
+                else:
+                    new_outputs.append(output)
+            new_outputs = tuple(new_outputs)
+        else:
+            new_outputs = outputs
+
+        return new_outputs
+
+    return wrapper

setup.py

@@ -45,7 +45,7 @@
 TESTS_REQUIRE = ['pytest', 'pytest-cov', 'pytest-mock', 'nose']
 EXTRAS_REQUIRE = {
     'optional': ['scipy', 'tables', 'numba', 'siphon', 'netcdf4',
-                 'ephem', 'cython'],
+                 'ephem', 'cython', 'pvfactors==0.1.5'],
     'doc': ['sphinx', 'ipython', 'sphinx_rtd_theme', 'numpydoc',
             'matplotlib'],
     'test': TESTS_REQUIRE