Implement reverse transposition using Perez-Driesse forward transposition

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.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`)
@@ -27,8 +30,10 @@ Documentation
 ~~~~~~~~~~~~~
 * Create :ref:`weatherdata` User's Guide page. (:pull:`1754`)
 * Fixed a plotting issue in the IV curve gallery example (:pull:`1895`)
+* Added two examples to demonstrate reverse transposition (:pull:`1907`)
 * Fixed `detect_clearsky` example in `clearsky.rst` (:issue:`1914`)
 
+
 Requirements
 ~~~~~~~~~~~~
 * Minimum version of scipy advanced from 1.4.0 to 1.5.0. (:issue:`1918`, :pull:`1919`)