Implement reverse transposition using Perez-Driesse forward transposition

docs/examples/irradiance-transposition/plot_rtranpose_limitations.py

@@ -0,0 +1,167 @@
+"""
+Explore the limitations of reverse transposition
+================================================
+
+Unfortunately, sometimes there is not a unique solution.
+
+Author: Anton Driesse
+
+"""
+
+# %%
+#
+# Introduction
+# ------------
+# 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:func:`pvlib.irradiance.rtranspose_driesse_2023` to estimate
+# the original GHI from POA global.
+#
+# References
+# ----------
+# .. [1] A. Driesse, A. Jensen, R. Perez, A Continuous Form of the Perez
+#     Diffuse Sky Model for Forward and Reverse Transposition, accepted
+#     for publication in the Solar Energy Journal.
+#
+
+import numpy as np
+
+import matplotlib
+import matplotlib.pyplot as plt
+
+from pvlib.irradiance import (erbs_driesse,
+                              get_total_irradiance,
+                              rtranspose_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 = rtranspose_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 = rtranspose_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,134 @@
+"""
+Reverse transposition using one year of hourly data
+===================================================
+
+With a brief look at accuracy and speed.
+
+Author: Anton Driesse
+
+"""
+
+# %%
+#
+# Introduction
+# ------------
+# In this example we start with a TMY file and calculate POA global irradiance.
+# Then, we use :py:func:`pvlib.irradiance.rtranspose_driesse_2023` to estimate
+# the original GHI from POA global.  Details of the method found in [1]_.
+#
+# References
+# ----------
+# .. [1] A. Driesse, A. Jensen, R. Perez, A Continuous Form of the Perez
+#     Diffuse Sky Model for Forward and Reverse Transposition, accepted
+#     for publication in the Solar Energy Journal.
+#
+
+import os
+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,
+                              rtranspose_driesse_2023,
+                              aoi,
+                              )
+
+from timeit import timeit
+
+# %%
+#
+# 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.
+#
+
+df['ghi_rev'] = rtranspose_driesse_2023(TILT, ORIENT,
+                                        solpos.apparent_zenith, solpos.azimuth,
+                                        df.poa_global,
+                                        dni_extra=df.dni_extra)
+
+# %%
+#
+# Let's see how long this takes...
+#
+
+
+def run_rtranspose():
+    rtranspose_driesse_2023(TILT, ORIENT,
+                            solpos.apparent_zenith, solpos.azimuth,
+                            df.poa_global,
+                            dni_extra=df.dni_extra)
+
+
+elapsed = timeit('run_rtranspose()', number=1, globals=globals())
+
+print('Elapsed time for reverse transposition: %.1f s' % elapsed)
+
+# %%
+#
+# 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.
+#
+
+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.rtranspose_driesse_2023

docs/sphinx/source/whatsnew/v0.10.3.rst

@@ -9,6 +9,9 @@ Enhancements
 ~~~~~~~~~~~~
 * Added the continuous Perez-Driesse transposition model.
   :py:func:`pvlib.irradiance.perez_driesse` (:issue:`1841`, :pull:`1876`)
+* Added a reverse transposition algorithm using the Perez-Driesse model.
+  :py:func:`pvlib.irradiance.rtranspose_driesse_2023`
+  (:issue:`1901`, :pull:`1907`)
 * :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance` and
   :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance_poa` now include
   shaded fraction in returned variables. (:pull:`1871`)
@@ -25,6 +28,7 @@ Testing
 Documentation
 ~~~~~~~~~~~~~
 * Fixed a plotting issue in the IV curve gallery example (:pull:`1895`)
+* Added two examples to demonstrate reverse transposition (:pull:`1907`)
 
 Contributors
 ~~~~~~~~~~~~