diff --git a/docs/sphinx/source/reference/irradiance/transposition.rst b/docs/sphinx/source/reference/irradiance/transposition.rst index 22136f0c58..d77e0d2801 100644 --- a/docs/sphinx/source/reference/irradiance/transposition.rst +++ b/docs/sphinx/source/reference/irradiance/transposition.rst @@ -15,4 +15,5 @@ Transposition models irradiance.klucher irradiance.reindl irradiance.king + irradiance.muneer irradiance.ghi_from_poa_driesse_2023 diff --git a/docs/sphinx/source/whatsnew/v0.11.1.rst b/docs/sphinx/source/whatsnew/v0.11.1.rst index db74f1362b..d6b414d018 100644 --- a/docs/sphinx/source/whatsnew/v0.11.1.rst +++ b/docs/sphinx/source/whatsnew/v0.11.1.rst @@ -30,6 +30,10 @@ Enhancements * Added function for calculating wind speed at different heights, :py:func:`pvlib.atmosphere.windspeed_powerlaw`. (:issue:`2118`, :pull:`2124`) +* Added function to determine sky diffuse irradiance on a tilted surface + using the Muneer model, + :py:func:`pvlib.irradiance.muneer`. + (:issue:`2117`, :pull:`2184`) Bug fixes ~~~~~~~~~ @@ -71,4 +75,3 @@ Contributors * Jose Meza (:ghuser:`JoseMezaMendieta`) * Bernat Nicolau (:ghuser:`BernatNicolau`) * Eduardo Sarquis (:ghuser:`EduardoSarquis`) - diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index 3b9c5cd0c0..369cb260f2 100644 --- a/pvlib/irradiance.py +++ b/pvlib/irradiance.py @@ -263,7 +263,7 @@ def get_total_irradiance(surface_tilt, surface_azimuth, dni, ghi, dhi, dni_extra=None, airmass=None, albedo=0.25, surface_type=None, model='isotropic', - model_perez='allsitescomposite1990'): + model_perez='allsitescomposite1990', b=5.73): r""" Determine total in-plane irradiance and its beam, sky diffuse and ground reflected components, using the specified sky diffuse irradiance model. @@ -280,6 +280,7 @@ def get_total_irradiance(surface_tilt, surface_azimuth, * king * perez * perez-driesse + * muneer Parameters ---------- @@ -309,9 +310,12 @@ def get_total_irradiance(surface_tilt, surface_azimuth, model : str, default 'isotropic' Irradiance model. Can be one of ``'isotropic'``, ``'klucher'``, ``'haydavies'``, ``'reindl'``, ``'king'``, ``'perez'``, - ``'perez-driesse'``. + ``'perez-driesse'`` and ``'muneer model_perez : str, default 'allsitescomposite1990' Used only if ``model='perez'``. See :py:func:`~pvlib.irradiance.perez`. + b : numeric, default 5.73 + Used only if ``model='muneer'``. + See :py:func:`~pvlib.irradiance.muneer`. Returns ------- @@ -321,19 +325,19 @@ def get_total_irradiance(surface_tilt, surface_azimuth, Notes ----- - Models ``'haydavies'``, ``'reindl'``, ``'perez'`` and ``'perez-driesse'`` - require ``'dni_extra'``. Values can be calculated using - :py:func:`~pvlib.irradiance.get_extra_radiation`. + Models ``'haydavies'``, ``'reindl'``, ``'perez'``, ``'perez-driesse'``, + and ``'muneer'`` require ``'dni_extra'``. Values of ``'dni_extra'`` + can be calculated using :py:func:`~pvlib.irradiance.get_extra_radiation`. The ``'perez'`` and ``'perez-driesse'`` models require relative airmass - (``airmass``) as input. If ``airmass`` is not provided, it is calculated + (``airmass``) as an input. If ``airmass`` is not provided, it is calculated using the defaults in :py:func:`~pvlib.atmosphere.get_relative_airmass`. """ poa_sky_diffuse = get_sky_diffuse( surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, dni, ghi, dhi, dni_extra=dni_extra, airmass=airmass, model=model, - model_perez=model_perez) + model_perez=model_perez, b=b) poa_ground_diffuse = get_ground_diffuse(surface_tilt, ghi, albedo, surface_type) @@ -346,7 +350,7 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth, dni, ghi, dhi, dni_extra=None, airmass=None, model='isotropic', - model_perez='allsitescomposite1990'): + model_perez='allsitescomposite1990', b=5.73): r""" Determine in-plane sky diffuse irradiance component using the specified sky diffuse irradiance model. @@ -359,6 +363,7 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, * king * perez * perez-driesse + * muneer Parameters ---------- @@ -383,9 +388,12 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, model : str, default 'isotropic' Irradiance model. Can be one of ``'isotropic'``, ``'klucher'``, ``'haydavies'``, ``'reindl'``, ``'king'``, ``'perez'``, - ``'perez-driesse'``. + ``'perez-driesse'``, ``'muneer'``. model_perez : str, default 'allsitescomposite1990' Used only if ``model='perez'``. See :py:func:`~pvlib.irradiance.perez`. + b : numeric, default 5.73 + Used only if ``model='muneer'``. + See :py:func:`~pvlib.irradiance.muneer`. Returns ------- @@ -395,14 +403,15 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, Raises ------ ValueError - If model is one of ``'haydavies'``, ``'reindl'``, ``'perez'``, or - ``'perez_driesse'`` and ``dni_extra`` is not specified. + If model is one of ``'haydavies'``, ``'reindl'``, ``'perez'``, + ``'perez_driesse'``, or ``'muneer'`` and ``dni_extra`` is not + specified. Notes ----- - Models ``'haydavies'``, ``'reindl'``, ``'perez'`` and ``'perez-driesse'`` - require ``'dni_extra'``. Values can be calculated using - :py:func:`~pvlib.irradiance.get_extra_radiation`. + Models ``'haydavies'``, ``'reindl'``, ``'perez'``, ``'perez-driesse'`` + and ``'muneer'`` require ``'dni_extra'``. Values can be calculated + using :py:func:`~pvlib.irradiance.get_extra_radiation`. The ``'Perez'`` transposition model features discontinuities in the predicted tilted diffuse irradiance due to relying on discrete input @@ -443,6 +452,9 @@ def get_sky_diffuse(surface_tilt, surface_azimuth, # perez_driesse will calculate its own airmass if needed sky = perez_driesse(surface_tilt, surface_azimuth, dhi, dni, dni_extra, solar_zenith, solar_azimuth, airmass) + elif model == 'muneer': + sky = muneer2004(surface_tilt, surface_azimuth, dhi, ghi, dni_extra, + solar_zenith, solar_azimuth) else: raise ValueError(f'invalid model selection {model}') @@ -994,6 +1006,252 @@ def king(surface_tilt, dhi, ghi, solar_zenith): return sky_diffuse +def muneer1990(surface_tilt, surface_azimuth, dhi, ghi, dni_extra, b=5.73, + solar_zenith=None, solar_azimuth=None, projection_ratio=None): + ''' + Determine sky diffuse irradiance on a tilted surface using the + Muneer model. + + This Muneer transposition model is described in [1]_. + + Parameters + ---------- + surface_tilt : numeric + Surface tilt angles in decimal degrees. ``surface_tilt`` must + be >=0 and <=180. The tilt angle is defined as degrees from horizontal + (e.g. surface facing up = 0, surface facing horizon = 90) + + surface_azimuth : numeric + Surface azimuth angles in decimal degrees. ``surface_azimuth`` + must be >=0 and <=360. The azimuth convention is defined as degrees + east of north (e.g. North = 0, South=180 East = 90, West = 270). + + dhi : numeric + Diffuse horizontal irradiance. [W/m^2] + + ghi : numeric + Global horizontal irradiance. [W/m^2] + + dni_extra : numeric + Extraterrestrial normal irradiance. [W/m^2] + + b : numeric, default 5.73 + Radiance distribution index, introduced by Moon and Spencer [2]_ + to model luminance distribution of overcast sky. [unitless] + Recommended values from [1]_: + + - isotropic: b = 0 + - shaded surface: b = 5.73 (default) + - sunlit surface, overcast sky: b = 1.68 + - sunlit surface, non-overcast sky: b = -0.62 + + solar_zenith : numeric + Solar apparent (refraction-corrected) zenith angles in decimal + degrees. Must supply ``solar_zenith`` and ``solar_azimuth`` or + supply ``projection_ratio``. + + solar_azimuth : numeric + Solar azimuth angles in decimal degrees. Must supply + ``solar_zenith`` and ``solar_azimuth`` or supply + ``projection_ratio``. + + projection_ratio : numeric, optional + Ratio of cosine of angle of incidence to cosine of solar zenith + angle. [unitless] + + Returns + ------- + poa_sky_diffuse : numeric + In-plane diffuse irradiance from the sky. [W/m^2] + + References + ---------- + .. [1] Muneer, T., 1990. Solar radiation model for Europe. Building + Services Engineering Research and Technology 11, 153-163. + :doi:`10.1177/014362449001100405` + + .. [2] Moon, P., and Spencer, D. E., 1942. Illumination from a + non-uniform sky. + Trans. Illum. Eng. Soc. (London) 37, 707-725. + ''' + + cos_solar_zenith = tools.cosd(solar_zenith) + # if necessary, calculate ratio of titled and horizontal beam irradiance + if projection_ratio is None: + cos_tt = aoi_projection(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth) + cos_tt = np.maximum(cos_tt, 0) # GH 526 + Rb = cos_tt / np.maximum(cos_solar_zenith, 0.01745) # GH 432 + else: + Rb = projection_ratio + + T_term1 = (1 + tools.cosd(surface_tilt)) * 0.5 + T_term2 = 2 * b / (np.pi * (3 + 2 * b)) + T_term3 = ( + tools.sind(surface_tilt) + - np.radians(surface_tilt) * tools.cosd(surface_tilt) + - np.pi * (1 - tools.cosd(surface_tilt)) * 0.5 + ) + T = T_term1 + T_term2 * T_term3 + + F = pvlib.irradiance.clearness_index( + ghi - dhi, solar_zenith, dni_extra) + + solar_elevation = np.pi/2 - np.radians(solar_zenith) + + numer_low = tools.sind(surface_tilt) * \ + tools.cosd(surface_azimuth - solar_azimuth) + denom_low = 0.1 - 0.008 * solar_elevation + + sky_diffuse_low = dhi*(T*(1-F) + F*(numer_low/denom_low)) + sky_diffuse_high = dhi*(T*(1-F) + F*Rb) + + low_elevation_condition = solar_elevation < 0.1 + sky_diffuse = np.where(low_elevation_condition, + sky_diffuse_low, sky_diffuse_high) + if isinstance(sky_diffuse_low, pd.Series): + sky_diffuse = pd.Series(sky_diffuse, + index=sky_diffuse_low.index) + + return sky_diffuse + + +def muneer2004(surface_tilt, surface_azimuth, dhi, ghi, dni_extra, + solar_zenith=None, solar_azimuth=None, projection_ratio=None, + location="southern_europe"): + ''' + Determine sky diffuse irradiance on a tilted surface using the + Muneer model. + + This Muneer transposition model is originally described in [1]_, and + improved for low solar elevations in [2]_. + + Parameters + ---------- + surface_tilt : numeric + Surface tilt angles in decimal degrees. ``surface_tilt`` must + be >=0 and <=180. The tilt angle is defined as degrees from horizontal + (e.g. surface facing up = 0, surface facing horizon = 90) + + surface_azimuth : numeric + Surface azimuth angles in decimal degrees. ``surface_azimuth`` + must be >=0 and <=360. The azimuth convention is defined as degrees + east of north (e.g. North = 0, South=180 East = 90, West = 270). + + dhi : numeric + Diffuse horizontal irradiance. [W/m^2] + + ghi : numeric + Global horizontal irradiance. [W/m^2] + + dni_extra : numeric + Extraterrestrial normal irradiance. [W/m^2] + + location : string, default 'northern_europe' + A string which selects the desired equation for calculating the second + term of T of Muneer model based on location. If location is not + provided as an input, the default, 'northern_europe' will be used. + All possible location selections are: + + * 'northern_europe' + * 'southern_europe' + * 'japan' + * 'globe' + + solar_zenith : numeric + Solar apparent (refraction-corrected) zenith angles in decimal + degrees. Must supply ``solar_zenith`` and ``solar_azimuth`` or + supply ``projection_ratio``. + + solar_azimuth : numeric + Solar azimuth angles in decimal degrees. Must supply + ``solar_zenith`` and ``solar_azimuth`` or supply + ``projection_ratio``. + + projection_ratio : numeric, optional + Ratio of cosine of angle of incidence to cosine of solar zenith + angle. [unitless] + + Returns + ------- + poa_sky_diffuse : numeric + In-plane diffuse irradiance from the sky. [W/m^2] + + References + ---------- + .. [1] Muneer, T., Gueymard, C., & Kambezidis, H., 2004. Hourly horizontal + irradiation and illuminance. 157-158 + :doi:`10.1016/B978-075065974-1/50011-5` + + .. [2] Gracia, A., & Huld, T., 2013. Performance comparison of different + models for the estimation of global irradiance on inclined surfaces: + validation of the model implemented in PVGIS. Centro Común de + Investigación, Instituto de Energía y Transporte. + :doi:`10.2790/91554` + + ''' + available_locations = ["northern_europe", "southern_europe", "japan", + "globe"] + if location not in available_locations: + raise ValueError + desired_index = available_locations.index(location) + + cos_solar_zenith = tools.cosd(solar_zenith) + # if necessary, calculate ratio of titled and horizontal beam irradiance + if projection_ratio is None: + cos_tt = aoi_projection(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth) + cos_tt = np.maximum(cos_tt, 0) # GH 526 + Rb = cos_tt / np.maximum(cos_solar_zenith, 0.01745) # GH 432 + else: + Rb = projection_ratio + + F = pvlib.irradiance.clearness_index( + ghi - dhi, solar_zenith, dni_extra) + + solar_elevation = np.pi/2 - np.radians(solar_zenith) + low_elevation_condition = solar_elevation < 0.1 + + aoi_ = aoi(surface_tilt, surface_azimuth, + solar_zenith, solar_azimuth) + overcast_condition = aoi_ >= 90 + F = np.where(overcast_condition, 0, F) + + T_term1 = (1 + tools.cosd(surface_tilt)) * 0.5 + + # term 2 depending on the selected location + T_term2_model0 = 0.00333 - 0.415 * F - 0.6987 * F ** 2 + T_term2_model1 = 0.00263 - 0.712 * F - 0.6883 * F ** 2 + T_term2_model2 = 0.08000 - 1.050 * F - 2.8400 * F ** 2 + T_term2_model3 = 0.04000 - 0.820 * F - 2.0260 * F ** 2 + T_term2 = np.array([T_term2_model0, T_term2_model1, + T_term2_model2, T_term2_model3]) + T_term2 = T_term2[desired_index] + T_term2 = np.where(overcast_condition, 0.25227, T_term2) + + T_term3 = ( + tools.sind(surface_tilt) + - np.radians(surface_tilt) * tools.cosd(surface_tilt) + - np.pi * (1 - tools.cosd(surface_tilt)) * 0.5 + ) + + T = T_term1 + T_term2 * T_term3 + + numer_low = tools.sind(surface_tilt) * \ + tools.cosd(surface_azimuth - solar_azimuth) + denom_low = 0.1 - 0.008 * solar_elevation + + sky_diffuse_low = dhi*(T*(1-F) + F*(numer_low/denom_low)) + sky_diffuse_high = dhi*(T*(1-F) + F*Rb) + + sky_diffuse = np.where(low_elevation_condition, + sky_diffuse_low, sky_diffuse_high) + if isinstance(ghi, pd.Series): + sky_diffuse = pd.Series(sky_diffuse, + index=ghi.index) + return sky_diffuse + + def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra, solar_zenith, solar_azimuth, airmass, model='allsitescomposite1990', return_components=False): diff --git a/pvlib/tests/test_irradiance.py b/pvlib/tests/test_irradiance.py index 0eb951f9a1..aa62f7b21e 100644 --- a/pvlib/tests/test_irradiance.py +++ b/pvlib/tests/test_irradiance.py @@ -286,6 +286,56 @@ def test_perez_driesse(irrad_data, ephem_data, dni_et, relative_airmass): assert_series_equal(out, expected, check_less_precise=2) +def test_muneer(irrad_data, ephem_data, dni_et): + projection_ratio = pd.Series(np.array( + [0, 0, 0.86399, 0.169725]), index=irrad_data.index) + dni_et = pd.Series(dni_et, index=irrad_data.index) + # pd.Series + out = irradiance.muneer(40, 180, irrad_data['dhi'], irrad_data['ghi'], + dni_et, solar_zenith=ephem_data['zenith'], + solar_azimuth=ephem_data['azimuth']) + out_Rb = irradiance.muneer(40, 180, irrad_data['dhi'], irrad_data['ghi'], + dni_et, + solar_zenith=ephem_data['zenith'], + solar_azimuth=ephem_data['azimuth'], + projection_ratio=projection_ratio) + # np.array + out_np = irradiance.muneer(40, 180, irrad_data['dhi'].values, + irrad_data['ghi'].values, dni_et.values, + solar_zenith=ephem_data['zenith'].values, + solar_azimuth=ephem_data['azimuth'].values) + out_np_Rb = irradiance.muneer(40, 180, irrad_data['dhi'].values, + irrad_data['ghi'].values, dni_et.values, + solar_zenith=ephem_data['zenith'].values, + solar_azimuth=ephem_data['azimuth'].values, + projection_ratio=projection_ratio.values) + # float + dhi_f = irrad_data['dhi'].values[2] + ghi_f = irrad_data['ghi'].values[2] + dni_et_f = dni_et.values[2] + solar_zenith_f = ephem_data['zenith'].values[2] + solar_azimuth_f = ephem_data['azimuth'].values[2] + projection_ratio_f = projection_ratio.values[2] + out_f = irradiance.muneer(40, 180, dhi_f, ghi_f, dni_et_f, + solar_zenith=solar_zenith_f, + solar_azimuth=solar_azimuth_f) + out_f_Rb = irradiance.muneer(40, 180, dhi_f, ghi_f, dni_et_f, + solar_zenith=solar_zenith_f, + solar_azimuth=solar_azimuth_f, + projection_ratio=projection_ratio_f) + + expected = pd.Series(np.array( + [0., 25.036, 100.759, 31.007]), + index=irrad_data.index) + expected_f = 100.759 + assert_series_equal(out, expected, check_less_precise=2) + assert_series_equal(out_Rb, expected, check_less_precise=2) + assert_almost_equal(out_np, expected.values, decimal=2) + assert_almost_equal(out_np_Rb, expected.values, decimal=2) + assert_almost_equal(out_f, expected_f, decimal=2) + assert_almost_equal(out_f_Rb, expected_f, decimal=2) + + def test_perez_driesse_airmass(irrad_data, ephem_data, dni_et): dni = irrad_data['dni'].copy() dni.iloc[2] = np.nan @@ -443,7 +493,8 @@ def test_perez_driesse_scalar(): @pytest.mark.parametrize('model', ['isotropic', 'klucher', 'haydavies', - 'reindl', 'king', 'perez', 'perez-driesse']) + 'reindl', 'king', 'perez', 'perez-driesse', + 'muneer']) def test_sky_diffuse_zenith_close_to_90(model): # GH 432 sky_diffuse = irradiance.get_sky_diffuse( @@ -495,7 +546,8 @@ def test_campbell_norman(): def test_get_total_irradiance(irrad_data, ephem_data, dni_et, relative_airmass): models = ['isotropic', 'klucher', - 'haydavies', 'reindl', 'king', 'perez', 'perez-driesse'] + 'haydavies', 'reindl', 'king', 'perez', 'perez-driesse', + 'muneer'] for model in models: total = irradiance.get_total_irradiance( @@ -514,7 +566,8 @@ def test_get_total_irradiance(irrad_data, ephem_data, dni_et, @pytest.mark.parametrize('model', ['isotropic', 'klucher', 'haydavies', 'reindl', 'king', - 'perez', 'perez-driesse']) + 'perez', 'perez-driesse', + 'muneer']) def test_get_total_irradiance_albedo( irrad_data, ephem_data, dni_et, relative_airmass, model): albedo = pd.Series(0.2, index=ephem_data.index) @@ -534,7 +587,8 @@ def test_get_total_irradiance_albedo( @pytest.mark.parametrize('model', ['isotropic', 'klucher', 'haydavies', 'reindl', 'king', - 'perez', 'perez-driesse']) + 'perez', 'perez-driesse', + 'muneer']) def test_get_total_irradiance_scalars(model): total = irradiance.get_total_irradiance( 32, 180,