Merge branch 'pvlib:main' into main

.github/workflows/pytest-remote-data.yml

@@ -56,7 +56,7 @@ jobs:
     strategy:
       fail-fast: false  # don't cancel other matrix jobs when one fails
       matrix:
-        python-version: [3.7, 3.8, 3.9, "3.10", "3.11"]
+        python-version: [3.7, 3.8, 3.9, "3.10", "3.11", "3.12"]
         suffix: ['']  # the alternative to "-min"
         include:
           - python-version: 3.7

.github/workflows/pytest.yml

@@ -12,7 +12,7 @@ jobs:
       fail-fast: false  # don't cancel other matrix jobs when one fails
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: [3.7, 3.8, 3.9, "3.10", "3.11"]
+        python-version: [3.7, 3.8, 3.9, "3.10", "3.11", "3.12"]
         environment-type: [conda, bare]
         suffix: ['']  # placeholder as an alternative to "-min"
         include:

ci/requirements-py3.12.yml

@@ -0,0 +1,28 @@
+name: test_env
+channels:
+    - defaults
+    - conda-forge
+dependencies:
+    - coveralls
+    - cython
+    - ephem
+    - h5py
+    # - numba  # not available for 3.12 as of 2023-12-12
+    - numpy >= 1.16.0
+    - pandas >= 0.25.0
+    - pip
+    - pytest
+    - pytest-cov
+    - pytest-mock
+    - requests-mock
+    - pytest-timeout
+    - pytest-rerunfailures
+    - pytest-remotedata
+    - python=3.12
+    - pytz
+    - requests
+    - scipy >= 1.5.0
+    - statsmodels
+    # - pip:
+        # - nrel-pysam>=2.0  # not available for 3.12 as of 2023-12-12
+        # - solarfactors  # required shapely<2 isn't available for 3.12

docs/examples/irradiance-transposition/plot_rtranpose_limitations.py

@@ -0,0 +1,181 @@
+"""
+Reverse transposition limitations
+====================================
+
+Unfortunately, sometimes there is not a unique solution.
+
+Author: Anton Driesse
+
+"""
+
+# %%
+#
+# Introduction
+# ------------
+# When irradiance is measured on a tilted plane, it is useful to be able to
+# estimate the GHI that produces the POA irradiance.
+# The estimation requires inverting a GHI-to-POA irradiance model,
+# which involves two parts:
+# a decomposition of GHI into direct and diffuse components,
+# and a transposition model that calculates the direct and diffuse irradiance
+# on the tilted plane.
+# Recovering GHI from POA irradiance is termed "reverse transposition."
+#
+# Unfortunately, for a given POA irradiance value, sometimes there is not a
+# unique solution for GHI.
+# Different GHI values can produce different combinations of direct and
+# diffuse irradiance that sum to the same POA irradiance value.
+#
+# In this example we look at a single point in time and consider a full range
+# of possible GHI and POA global values as shown in figures 3 and 4 of [1]_.
+# Then we use :py:meth:`pvlib.irradiance.ghi_from_poa_driesse_2023` to estimate
+# the original GHI from POA global.
+#
+# References
+# ----------
+# .. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the
+#     Perez diffuse sky model for forward and reverse transposition.
+#     Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093`
+#
+
+import numpy as np
+
+import matplotlib
+import matplotlib.pyplot as plt
+
+from pvlib.irradiance import (erbs_driesse,
+                              get_total_irradiance,
+                              ghi_from_poa_driesse_2023,
+                              )
+
+matplotlib.rcParams['axes.grid'] = True
+
+# %%
+#
+# Define the conditions that were used for figure 3 in [1]_.
+#
+
+dni_extra = 1366.1
+albedo = 0.25
+surface_tilt = 40
+surface_azimuth = 180
+
+solar_azimuth = 82
+solar_zenith = 75
+
+# %%
+#
+# Define a range of possible GHI values and calculate the corresponding
+# POA global.  First estimate DNI and DHI using the Erbs-Driesse model, then
+# transpose using the Perez-Driesse model.
+#
+
+ghi = np.linspace(0, 500, 100+1)
+
+erbsout = erbs_driesse(ghi, solar_zenith, dni_extra=dni_extra)
+
+dni = erbsout['dni']
+dhi = erbsout['dhi']
+
+irrads = get_total_irradiance(surface_tilt, surface_azimuth,
+                              solar_zenith, solar_azimuth,
+                              dni, ghi, dhi,
+                              dni_extra,
+                              model='perez-driesse')
+
+poa_global = irrads['poa_global']
+
+# %%
+#
+# Suppose you measure that POA global is 200 W/m2. What would GHI be?
+#
+
+poa_test = 200
+
+ghi_hat = ghi_from_poa_driesse_2023(surface_tilt, surface_azimuth,
+                                    solar_zenith, solar_azimuth,
+                                    poa_test,
+                                    dni_extra,
+                                    full_output=False)
+
+print('Estimated GHI: %.2f W/m².' % ghi_hat)
+
+# %%
+#
+# Show this result on the graph of all possible combinations of GHI and POA.
+#
+
+plt.figure()
+plt.plot(ghi, poa_global, 'k-')
+plt.axvline(ghi_hat, color='g', lw=1)
+plt.axhline(poa_test, color='g', lw=1)
+plt.plot(ghi_hat, poa_test, 'gs')
+plt.annotate('GHI=%.2f' % (ghi_hat),
+             xy=(ghi_hat-2, 200+2),
+             xytext=(ghi_hat-20, 200+20),
+             ha='right',
+             arrowprops={'arrowstyle': 'simple'})
+plt.xlim(0, 500)
+plt.ylim(0, 250)
+plt.xlabel('GHI [W/m²]')
+plt.ylabel('POA [W/m²]')
+plt.show()
+
+# %%
+#
+# Now change the solar azimuth to match the conditions for figure 4 in [1]_.
+#
+
+solar_azimuth = 76
+
+# %%
+#
+# Again, estimate DNI and DHI using the Erbs-Driesse model, then
+# transpose using the Perez-Driesse model.
+#
+
+erbsout = erbs_driesse(ghi, solar_zenith, dni_extra=dni_extra)
+
+dni = erbsout['dni']
+dhi = erbsout['dhi']
+
+irrads = get_total_irradiance(surface_tilt, surface_azimuth,
+                              solar_zenith, solar_azimuth,
+                              dni, ghi, dhi,
+                              dni_extra,
+                              model='perez-driesse')
+
+poa_global = irrads['poa_global']
+
+# %%
+#
+# Now reverse transpose all the POA values and observe that the original
+# GHI cannot always be found.  There is a range of POA values that
+# maps to three possible GHI values, and there is not enough information
+# to choose one of them.  Sometimes we get lucky and the right one comes
+# out, other times not.
+#
+
+result = ghi_from_poa_driesse_2023(surface_tilt, surface_azimuth,
+                                   solar_zenith, solar_azimuth,
+                                   poa_global,
+                                   dni_extra,
+                                   full_output=True,
+                                   )
+
+ghi_hat, conv, niter = result
+correct = np.isclose(ghi, ghi_hat, atol=0.01)
+
+plt.figure()
+plt.plot(np.where(correct, ghi, np.nan), np.where(correct, poa_global, np.nan),
+         'g.', label='correct GHI found')
+plt.plot(ghi[~correct], poa_global[~correct], 'r.', label='unreachable GHI')
+plt.plot(ghi[~conv], poa_global[~conv], 'm.', label='out of range (kt > 1.25)')
+plt.axhspan(88, 103, color='y', alpha=0.25, label='problem region')
+
+plt.xlim(0, 500)
+plt.ylim(0, 250)
+plt.xlabel('GHI [W/m²]')
+plt.ylabel('POA [W/m²]')
+plt.legend()
+plt.show()

docs/examples/irradiance-transposition/plot_rtranpose_year.py

@@ -0,0 +1,152 @@
+"""
+Reverse transposition using one year of hourly data
+===================================================
+
+With a brief look at accuracy and speed.
+
+Author: Anton Driesse
+
+"""
+# %%
+#
+# Introduction
+# ------------
+# When irradiance is measured on a tilted plane, it is useful to be able to
+# estimate the GHI that produces the POA irradiance.
+# The estimation requires inverting a GHI-to-POA irradiance model,
+# which involves two parts:
+# a decomposition of GHI into direct and diffuse components,
+# and a transposition model that calculates the direct and diffuse
+# irradiance on the tilted plane.
+# Recovering GHI from POA irradiance is termed "reverse transposition."
+#
+# In this example we start with a TMY file and calculate POA global irradiance.
+# Then we use :py:meth:`pvlib.irradiance.ghi_from_poa_driesse_2023` to estimate
+# the original GHI from POA global.  Details of the method found in [1]_.
+#
+# Another method for reverse tranposition called GTI-DIRINT is also
+# available in pvlib python (:py:meth:`pvlib.irradiance.gti_dirint`).
+# More information is available in [2]_.
+#
+# References
+# ----------
+# .. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the
+#     Perez diffuse sky model for forward and reverse transposition.
+#     Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093`
+#
+# .. [2] B. Marion, A model for deriving the direct normal and
+#        diffuse horizontal irradiance from the global tilted
+#        irradiance, Solar Energy 122, 1037-1046.
+#        :doi:`10.1016/j.solener.2015.10.024`
+
+import os
+import time
+import pandas as pd
+
+import matplotlib.pyplot as plt
+
+import pvlib
+from pvlib import iotools, location
+from pvlib.irradiance import (get_extra_radiation,
+                              get_total_irradiance,
+                              ghi_from_poa_driesse_2023,
+                              aoi,
+                              )
+
+# %%
+#
+# Read a TMY3 file containing weather data and select needed columns.
+#
+
+PVLIB_DIR = pvlib.__path__[0]
+DATA_FILE = os.path.join(PVLIB_DIR, 'data', '723170TYA.CSV')
+
+tmy, metadata = iotools.read_tmy3(DATA_FILE, coerce_year=1990,
+                                  map_variables=True)
+
+df = pd.DataFrame({'ghi': tmy['ghi'], 'dhi': tmy['dhi'], 'dni': tmy['dni'],
+                   'temp_air': tmy['temp_air'],
+                   'wind_speed': tmy['wind_speed'],
+                   })
+
+# %%
+#
+# Shift the timestamps to the middle of the hour and calculate sun positions.
+#
+
+df.index = df.index - pd.Timedelta(minutes=30)
+
+loc = location.Location.from_tmy(metadata)
+solpos = loc.get_solarposition(df.index)
+
+# %%
+#
+# Estimate global irradiance on a fixed-tilt array (forward transposition).
+# The array is tilted 30 degrees and oriented 30 degrees east of south.
+#
+
+TILT = 30
+ORIENT = 150
+
+df['dni_extra'] = get_extra_radiation(df.index)
+
+total_irrad = get_total_irradiance(TILT, ORIENT,
+                                   solpos.apparent_zenith,
+                                   solpos.azimuth,
+                                   df.dni, df.ghi, df.dhi,
+                                   dni_extra=df.dni_extra,
+                                   model='perez-driesse')
+
+df['poa_global'] = total_irrad.poa_global
+df['aoi'] = aoi(TILT, ORIENT, solpos.apparent_zenith, solpos.azimuth)
+
+# %%
+#
+# Now estimate ghi from poa_global using reverse transposition.
+# The algorithm uses a simple bisection search, which is quite slow
+# because scipy doesn't offer a vectorized version (yet).
+# For this reason we'll process a random sample of 1000 timestamps
+# rather than the whole year.
+#
+
+df = df[df.ghi > 0].sample(n=1000)
+solpos = solpos.reindex(df.index)
+
+start = time.process_time()
+
+df['ghi_rev'] = ghi_from_poa_driesse_2023(TILT, ORIENT,
+                                          solpos.apparent_zenith,
+                                          solpos.azimuth,
+                                          df.poa_global,
+                                          dni_extra=df.dni_extra)
+finish = time.process_time()
+
+print('Elapsed time for reverse transposition: %.1f s' % (finish - start))
+
+# %%
+#
+# This graph shows the reverse transposed values vs. the original values.
+# The markers are color-coded by angle-of-incidence to show that
+# errors occur primarily with incidence angle approaching 90° and beyond.
+#
+# Note that the results look particularly good because the POA values
+# were calculated using the same models as used in reverse transposition.
+# This isn't cheating though.  It's a way of ensuring that the errors
+# we see are really due to the reverse transposition algorithm.
+# Expect to see larger errors with real-word POA measurements
+# because errors from forward and reverse transposition will both be present.
+#
+
+df = df.sort_values('aoi')
+
+plt.figure()
+plt.gca().grid(True, alpha=.5)
+pc = plt.scatter(df['ghi'], df['ghi_rev'], c=df['aoi'], s=15,
+                 cmap='jet', vmin=60, vmax=120)
+plt.colorbar(label='AOI [°]')
+pc.set_alpha(0.5)
+
+plt.xlabel('GHI original [W/m²]')
+plt.ylabel('GHI from POA [W/m²]')
+
+plt.show()

docs/sphinx/source/reference/irradiance/transposition.rst

@@ -15,3 +15,4 @@ Transposition models
    irradiance.klucher
    irradiance.reindl
    irradiance.king
+   irradiance.ghi_from_poa_driesse_2023