Skip to content

Commit c18a004

Browse files
adriessecwhanseAdamRJensen
authored
Implement reverse transposition using Perez-Driesse forward transposition (#1907)
* Add reverse transposition function and two helpers. * Add missing import. * Add to docs under transposition until a better place is found/made. * Minor doc string fixes. * Add full_output option similar to newton(). * First example for reverse transposition. * Second example for reverse transposition. * Some improvements to the two examples. * Placate flake8. * Add tests for reverse transpostion. * Refine examples and fix test. * Update whatsnew. * Refine examples. * Try to get rid of matplotlib warning in example. * Remove unused import. * Update pvlib/irradiance.py Co-authored-by: Cliff Hansen <[email protected]> * Improve examples based on reviews. * Settle conflict. * Try again. * Remove one space. * Final(?) changes. * Update reference in erbs_driesse(). * Fix links in examples. * Update docs/examples/irradiance-transposition/plot_rtranpose_limitations.py Co-authored-by: Adam R. Jensen <[email protected]> * Address review comments. * Final renames. --------- Co-authored-by: Cliff Hansen <[email protected]> Co-authored-by: Adam R. Jensen <[email protected]>
1 parent 300aedc commit c18a004

File tree

6 files changed

+554
-5
lines changed

6 files changed

+554
-5
lines changed

‎docs/examples/irradiance-transposition/plot_rtranpose_limitations.py

Lines changed: 181 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,181 @@
1+
"""
2+
Reverse transposition limitations
3+
====================================
4+
5+
Unfortunately, sometimes there is not a unique solution.
6+
7+
Author: Anton Driesse
8+
9+
"""
10+
11+
# %%
12+
#
13+
# Introduction
14+
# ------------
15+
# When irradiance is measured on a tilted plane, it is useful to be able to
16+
# estimate the GHI that produces the POA irradiance.
17+
# The estimation requires inverting a GHI-to-POA irradiance model,
18+
# which involves two parts:
19+
# a decomposition of GHI into direct and diffuse components,
20+
# and a transposition model that calculates the direct and diffuse irradiance
21+
# on the tilted plane.
22+
# Recovering GHI from POA irradiance is termed "reverse transposition."
23+
#
24+
# Unfortunately, for a given POA irradiance value, sometimes there is not a
25+
# unique solution for GHI.
26+
# Different GHI values can produce different combinations of direct and
27+
# diffuse irradiance that sum to the same POA irradiance value.
28+
#
29+
# In this example we look at a single point in time and consider a full range
30+
# of possible GHI and POA global values as shown in figures 3 and 4 of [1]_.
31+
# Then we use :py:meth:`pvlib.irradiance.ghi_from_poa_driesse_2023` to estimate
32+
# the original GHI from POA global.
33+
#
34+
# References
35+
# ----------
36+
# .. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the
37+
# Perez diffuse sky model for forward and reverse transposition.
38+
# Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093`
39+
#
40+
41+
import numpy as np
42+
43+
import matplotlib
44+
import matplotlib.pyplot as plt
45+
46+
from pvlib.irradiance import (erbs_driesse,
47+
get_total_irradiance,
48+
ghi_from_poa_driesse_2023,
49+
)
50+
51+
matplotlib.rcParams['axes.grid'] = True
52+
53+
# %%
54+
#
55+
# Define the conditions that were used for figure 3 in [1]_.
56+
#
57+
58+
dni_extra = 1366.1
59+
albedo = 0.25
60+
surface_tilt = 40
61+
surface_azimuth = 180
62+
63+
solar_azimuth = 82
64+
solar_zenith = 75
65+
66+
# %%
67+
#
68+
# Define a range of possible GHI values and calculate the corresponding
69+
# POA global. First estimate DNI and DHI using the Erbs-Driesse model, then
70+
# transpose using the Perez-Driesse model.
71+
#
72+
73+
ghi = np.linspace(0, 500, 100+1)
74+
75+
erbsout = erbs_driesse(ghi, solar_zenith, dni_extra=dni_extra)
76+
77+
dni = erbsout['dni']
78+
dhi = erbsout['dhi']
79+
80+
irrads = get_total_irradiance(surface_tilt, surface_azimuth,
81+
solar_zenith, solar_azimuth,
82+
dni, ghi, dhi,
83+
dni_extra,
84+
model='perez-driesse')
85+
86+
poa_global = irrads['poa_global']
87+
88+
# %%
89+
#
90+
# Suppose you measure that POA global is 200 W/m2. What would GHI be?
91+
#
92+
93+
poa_test = 200
94+
95+
ghi_hat = ghi_from_poa_driesse_2023(surface_tilt, surface_azimuth,
96+
solar_zenith, solar_azimuth,
97+
poa_test,
98+
dni_extra,
99+
full_output=False)
100+
101+
print('Estimated GHI: %.2f W/m².' % ghi_hat)
102+
103+
# %%
104+
#
105+
# Show this result on the graph of all possible combinations of GHI and POA.
106+
#
107+
108+
plt.figure()
109+
plt.plot(ghi, poa_global, 'k-')
110+
plt.axvline(ghi_hat, color='g', lw=1)
111+
plt.axhline(poa_test, color='g', lw=1)
112+
plt.plot(ghi_hat, poa_test, 'gs')
113+
plt.annotate('GHI=%.2f' % (ghi_hat),
114+
xy=(ghi_hat-2, 200+2),
115+
xytext=(ghi_hat-20, 200+20),
116+
ha='right',
117+
arrowprops={'arrowstyle': 'simple'})
118+
plt.xlim(0, 500)
119+
plt.ylim(0, 250)
120+
plt.xlabel('GHI [W/m²]')
121+
plt.ylabel('POA [W/m²]')
122+
plt.show()
123+
124+
# %%
125+
#
126+
# Now change the solar azimuth to match the conditions for figure 4 in [1]_.
127+
#
128+
129+
solar_azimuth = 76
130+
131+
# %%
132+
#
133+
# Again, estimate DNI and DHI using the Erbs-Driesse model, then
134+
# transpose using the Perez-Driesse model.
135+
#
136+
137+
erbsout = erbs_driesse(ghi, solar_zenith, dni_extra=dni_extra)
138+
139+
dni = erbsout['dni']
140+
dhi = erbsout['dhi']
141+
142+
irrads = get_total_irradiance(surface_tilt, surface_azimuth,
143+
solar_zenith, solar_azimuth,
144+
dni, ghi, dhi,
145+
dni_extra,
146+
model='perez-driesse')
147+
148+
poa_global = irrads['poa_global']
149+
150+
# %%
151+
#
152+
# Now reverse transpose all the POA values and observe that the original
153+
# GHI cannot always be found. There is a range of POA values that
154+
# maps to three possible GHI values, and there is not enough information
155+
# to choose one of them. Sometimes we get lucky and the right one comes
156+
# out, other times not.
157+
#
158+
159+
result = ghi_from_poa_driesse_2023(surface_tilt, surface_azimuth,
160+
solar_zenith, solar_azimuth,
161+
poa_global,
162+
dni_extra,
163+
full_output=True,
164+
)
165+
166+
ghi_hat, conv, niter = result
167+
correct = np.isclose(ghi, ghi_hat, atol=0.01)
168+
169+
plt.figure()
170+
plt.plot(np.where(correct, ghi, np.nan), np.where(correct, poa_global, np.nan),
171+
'g.', label='correct GHI found')
172+
plt.plot(ghi[~correct], poa_global[~correct], 'r.', label='unreachable GHI')
173+
plt.plot(ghi[~conv], poa_global[~conv], 'm.', label='out of range (kt > 1.25)')
174+
plt.axhspan(88, 103, color='y', alpha=0.25, label='problem region')
175+
176+
plt.xlim(0, 500)
177+
plt.ylim(0, 250)
178+
plt.xlabel('GHI [W/m²]')
179+
plt.ylabel('POA [W/m²]')
180+
plt.legend()
181+
plt.show()

‎docs/examples/irradiance-transposition/plot_rtranpose_year.py

Lines changed: 152 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,152 @@
1+
"""
2+
Reverse transposition using one year of hourly data
3+
===================================================
4+
5+
With a brief look at accuracy and speed.
6+
7+
Author: Anton Driesse
8+
9+
"""
10+
# %%
11+
#
12+
# Introduction
13+
# ------------
14+
# When irradiance is measured on a tilted plane, it is useful to be able to
15+
# estimate the GHI that produces the POA irradiance.
16+
# The estimation requires inverting a GHI-to-POA irradiance model,
17+
# which involves two parts:
18+
# a decomposition of GHI into direct and diffuse components,
19+
# and a transposition model that calculates the direct and diffuse
20+
# irradiance on the tilted plane.
21+
# Recovering GHI from POA irradiance is termed "reverse transposition."
22+
#
23+
# In this example we start with a TMY file and calculate POA global irradiance.
24+
# Then we use :py:meth:`pvlib.irradiance.ghi_from_poa_driesse_2023` to estimate
25+
# the original GHI from POA global. Details of the method found in [1]_.
26+
#
27+
# Another method for reverse tranposition called GTI-DIRINT is also
28+
# available in pvlib python (:py:meth:`pvlib.irradiance.gti_dirint`).
29+
# More information is available in [2]_.
30+
#
31+
# References
32+
# ----------
33+
# .. [1] Driesse, A., Jensen, A., Perez, R., 2024. A Continuous form of the
34+
# Perez diffuse sky model for forward and reverse transposition.
35+
# Solar Energy vol. 267. :doi:`10.1016/j.solener.2023.112093`
36+
#
37+
# .. [2] B. Marion, A model for deriving the direct normal and
38+
# diffuse horizontal irradiance from the global tilted
39+
# irradiance, Solar Energy 122, 1037-1046.
40+
# :doi:`10.1016/j.solener.2015.10.024`
41+
42+
import os
43+
import time
44+
import pandas as pd
45+
46+
import matplotlib.pyplot as plt
47+
48+
import pvlib
49+
from pvlib import iotools, location
50+
from pvlib.irradiance import (get_extra_radiation,
51+
get_total_irradiance,
52+
ghi_from_poa_driesse_2023,
53+
aoi,
54+
)
55+
56+
# %%
57+
#
58+
# Read a TMY3 file containing weather data and select needed columns.
59+
#
60+
61+
PVLIB_DIR = pvlib.__path__[0]
62+
DATA_FILE = os.path.join(PVLIB_DIR, 'data', '723170TYA.CSV')
63+
64+
tmy, metadata = iotools.read_tmy3(DATA_FILE, coerce_year=1990,
65+
map_variables=True)
66+
67+
df = pd.DataFrame({'ghi': tmy['ghi'], 'dhi': tmy['dhi'], 'dni': tmy['dni'],
68+
'temp_air': tmy['temp_air'],
69+
'wind_speed': tmy['wind_speed'],
70+
})
71+
72+
# %%
73+
#
74+
# Shift the timestamps to the middle of the hour and calculate sun positions.
75+
#
76+
77+
df.index = df.index - pd.Timedelta(minutes=30)
78+
79+
loc = location.Location.from_tmy(metadata)
80+
solpos = loc.get_solarposition(df.index)
81+
82+
# %%
83+
#
84+
# Estimate global irradiance on a fixed-tilt array (forward transposition).
85+
# The array is tilted 30 degrees and oriented 30 degrees east of south.
86+
#
87+
88+
TILT = 30
89+
ORIENT = 150
90+
91+
df['dni_extra'] = get_extra_radiation(df.index)
92+
93+
total_irrad = get_total_irradiance(TILT, ORIENT,
94+
solpos.apparent_zenith,
95+
solpos.azimuth,
96+
df.dni, df.ghi, df.dhi,
97+
dni_extra=df.dni_extra,
98+
model='perez-driesse')
99+
100+
df['poa_global'] = total_irrad.poa_global
101+
df['aoi'] = aoi(TILT, ORIENT, solpos.apparent_zenith, solpos.azimuth)
102+
103+
# %%
104+
#
105+
# Now estimate ghi from poa_global using reverse transposition.
106+
# The algorithm uses a simple bisection search, which is quite slow
107+
# because scipy doesn't offer a vectorized version (yet).
108+
# For this reason we'll process a random sample of 1000 timestamps
109+
# rather than the whole year.
110+
#
111+
112+
df = df[df.ghi > 0].sample(n=1000)
113+
solpos = solpos.reindex(df.index)
114+
115+
start = time.process_time()
116+
117+
df['ghi_rev'] = ghi_from_poa_driesse_2023(TILT, ORIENT,
118+
solpos.apparent_zenith,
119+
solpos.azimuth,
120+
df.poa_global,
121+
dni_extra=df.dni_extra)
122+
finish = time.process_time()
123+
124+
print('Elapsed time for reverse transposition: %.1f s' % (finish - start))
125+
126+
# %%
127+
#
128+
# This graph shows the reverse transposed values vs. the original values.
129+
# The markers are color-coded by angle-of-incidence to show that
130+
# errors occur primarily with incidence angle approaching 90° and beyond.
131+
#
132+
# Note that the results look particularly good because the POA values
133+
# were calculated using the same models as used in reverse transposition.
134+
# This isn't cheating though. It's a way of ensuring that the errors
135+
# we see are really due to the reverse transposition algorithm.
136+
# Expect to see larger errors with real-word POA measurements
137+
# because errors from forward and reverse transposition will both be present.
138+
#
139+
140+
df = df.sort_values('aoi')
141+
142+
plt.figure()
143+
plt.gca().grid(True, alpha=.5)
144+
pc = plt.scatter(df['ghi'], df['ghi_rev'], c=df['aoi'], s=15,
145+
cmap='jet', vmin=60, vmax=120)
146+
plt.colorbar(label='AOI [°]')
147+
pc.set_alpha(0.5)
148+
149+
plt.xlabel('GHI original [W/m²]')
150+
plt.ylabel('GHI from POA [W/m²]')
151+
152+
plt.show()

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,3 +15,4 @@ Transposition models
1515
irradiance.klucher
1616
irradiance.reindl
1717
irradiance.king
18+
irradiance.ghi_from_poa_driesse_2023

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

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,9 @@ Enhancements
99
~~~~~~~~~~~~
1010
* Added the continuous Perez-Driesse transposition model.
1111
:py:func:`pvlib.irradiance.perez_driesse` (:issue:`1841`, :pull:`1876`)
12+
* Added a reverse transposition algorithm using the Perez-Driesse model.
13+
:py:func:`pvlib.irradiance.ghi_from_poa_driesse_2023`
14+
(:issue:`1901`, :pull:`1907`)
1215
* :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance` and
1316
:py:func:`pvlib.bifacial.infinite_sheds.get_irradiance_poa` now include
1417
shaded fraction in returned variables. (:pull:`1871`)
@@ -27,8 +30,10 @@ Documentation
2730
~~~~~~~~~~~~~
2831
* Create :ref:`weatherdata` User's Guide page. (:pull:`1754`)
2932
* Fixed a plotting issue in the IV curve gallery example (:pull:`1895`)
33+
* Added two examples to demonstrate reverse transposition (:pull:`1907`)
3034
* Fixed `detect_clearsky` example in `clearsky.rst` (:issue:`1914`)
3135

36+
3237
Requirements
3338
~~~~~~~~~~~~
3439
* Minimum version of scipy advanced from 1.4.0 to 1.5.0. (:issue:`1918`, :pull:`1919`)

0 commit comments

Comments
 (0)