<jupyter_start><jupyter_text>Solar Decomposition ModelsKnowing the direct or beam normal irradiance (DNI) is useful for many solar and energy applications, e.g., calculating the yield of solar concetrating power systems or determining the irradiance on an inclined surface. However, often only global horizontal irradiance (GHI) is available, as it is much cheaper to measure than the individual components and most satelitte radiation models only calculate GHI. Many applications, therefore, require applying **solar radiation decomposition models** to obtain the direct beam irradiance. The majority of decomposition models work by applying an empirical correlations between GHI, diffuse horizontal irradiance (DHI) and DNI. In this section, we present the implementation in Python of several solar decomposition models. *** Content by <NAME> and <NAME>   The role of solar decomposition modelsSolar decomposition models are some of the most frequently used algorithms, and are often used during the first part of a larger analysis. As described earlier, decomposition models permit estimating the beam or direct normal irradiance (DNI) from GHI measurements, often using relationships between the clearness index ($k_t$) and the diffuse fraction of solar irradiance ($k_d$). The clearness index is defined as:$$ k_t = \frac{GHI}{E_a \cdot cos(\theta_Z)}$$where $E_a$ is the extraterrestrial irradiance, and $\theta_Z$ is the solar zenith angle. The diffuse fraction is defined as: $$ k_d = \frac{DHI}{GHI} $$The role of decomposition models can be graphically illustrated in the diagram below:```{image} ../notebooks/graphics/decomposition_models.png:alt: decomposition_flowchart:class: bg-white mb-1:width: 500px:align: center``` Multiple models have been proposed for the decomposition of global horizontal radiation/irradiance into its direct normal and diffuse horizontal components in the scientific literature. Below, the implementation in Python of some of these methods will be shown. The following methods will be modelled and compared:- DISC model (Maxwell, 1987)- DIRINT model (Pérez et al., 1992)- DIRINDEX model (Pérez et al., 2002)- Erbs model (Erbs et al., 1982)The four decomposition models, among some others, are implemented in the open-source library [pvlib python](https://pvlib-python.readthedocs.io/en/stable/api.htmldni-estimation-models).Before applying the models in practice, an example irradiance data set will be retrieved.  Data preparationTo illustrate the implementation of each of these models, we will use real irradiance observations. The decomposition models will be used to calculate DNI and DIh from measured GHI, and the results will then be compared to measurements of DNI and DHI. The first step is getting and preparing the measurement data to be used.<jupyter_code># Import all the necessary libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pvlib<jupyter_output><empty_output><jupyter_text>Retrieving weather dataThe irradiance data used is from the [BSRN](https://bsrn.awi.de/) Cabauw station operated by the Royal Netherlands Meteorological Institute (KNMI). The data is retrieved from the BSRN FTP server using the `get_bsrn` function from pvlib-python. The Cabauw station is located at latitude 51.9711° and longitude 4.9267°, at an elevation of 0 m AMSL. The reason for choosing the Cabauw station is that the data has already been extensively quality checked, hence in this example this step is skipped. The data used are 1-minute global, direct, and diffuse irradiance measurements between January and July 2021, which are resampled to hourly observations.<jupyter_code># Load the BSRN credentials from the Secret environment variables
# This cell is hidden in the displayed notebook
import os
#bsrn_username = os.environ["BSRN_FTP_USERNAME"]
#bsrn_password = os.environ["BSRN_FTP_PASSWORD"]
bsrn_username = 'bsrnftp'
bsrn_password = '<PASSWORD>'
sodapro_username = '<EMAIL>'
# Retrieving the measurement data
df, meta = pvlib.iotools.get_bsrn(
    station='CAB',
    start=pd.Timestamp('2021-01-01'),
    end=pd.Timestamp('2021-07-31'),
    username=bsrn_username,
    password=<PASSWORD>)

# Resample the 1-minute data to hourly
df = df.resample('1h').mean()

# Display the first five rows of the dataset
df.head()<jupyter_output><empty_output><jupyter_text>Besides GHI, DNI and DHI, this station has additional weather-related variables available (e.g., air temperature, humidity, and wind speed), as well as sensor-specific data. The details of the sensors and measurements available in each station can be found under the section "Instrument history and meta data" of each MIDC's station webpage. Let's list the variables available and then subset those that we will use in this example:<jupyter_code># Visualize variables avilable
df.columns<jupyter_output><empty_output><jupyter_text>It is often helpful to visualize the data in order to check that the measurements are within the expected ones. For example, a plot for a single day (18th July 2021):<jupyter_code>df.loc['2021-07-18', ['ghi','dhi','dni']].plot()
plt.title('Solar irradiance components, 18th July 2021')
plt.ylabel('Irradiance [W/m$^2$]')
plt.legend()<jupyter_output><empty_output><jupyter_text>For an overall representation of the data, it's useful to inspect the the key statistical metrics for the data:<jupyter_code>df[['ghi','dni','dhi', 'temp_air', 'relative_humidity', 'pressure']].describe()<jupyter_output><empty_output><jupyter_text>**Estimating other input variables:**The 3 solar decomposition models that we will implement in this section (i.e., DISC, DIRINT and Erbs models) have different inputs to be implemented. **Solar position data** is a common input for the 3 models. The DISC and DIRINT models accept as a possible input the **dew (wet-bulb) temperature** ($T_d$). This variable is not directly available from our weather observations; however, it can be easily estimated using a simple conversion that has relative humidity ($RH$) and the ambient temperature ($T_a$) as input. This simple approximation is given by the equation below (Lawrence, 2005): $$ T_{d} = T_a (^{\circ} C) - \frac{100-RH (\%)}{5} $$ In Python, this equation could be implemented as follows:<jupyter_code>df['temp_dew'] = df['temp_air']-((100-df['relative_humidity'])/5)<jupyter_output><empty_output><jupyter_text>The solar position for the entire period can be calculated using the `pvlib.solarposition.get_solarposition` function which supports a number of different solar position algorithms. The default algorithm is an implementation of the NREL SPA algorithm.<jupyter_code># Define a pvlib location object for the site
location = pvlib.location.Location(latitude=51.9711, longitude=4.9267, altitude=0, name='Cabauw')

# Calculate the solar position for each index in the datarame for the specific location
# The 30 minute shift is applied in order to calculate the solar position at the middle of the hour!
solpos = location.get_solarposition(df.index+pd.Timedelta(minutes=30))
solpos.index = solpos.index - pd.Timedelta(minutes=30)

# Let's have a look to the output
solpos.head()<jupyter_output><empty_output><jupyter_text>DISC modelThe DISC model (Maxwell, 1987) requires the following inputs: GHI, solar zenith angle, and the atmospheric pressure, which is used to estimate the absolute (pressure-corrected) airmass.<jupyter_code># DISC method estimated with absolute airmass as input
disc = pvlib.irradiance.disc(
    ghi=df['ghi'],
    solar_zenith=solpos['apparent_zenith'],
    datetime_or_doy=df.index,
    pressure=df['pressure']*100) # Site pressure in Pascal<jupyter_output><empty_output><jupyter_text>The functoin returns the DNI, the clearness index ($k_t$) and the airmass estimation:<jupyter_code>disc.columns<jupyter_output><empty_output><jupyter_text>DIRINT modelAn implementatin of the DIRINT model (Pérez et al., 1992)  is also available in pvlib-python. The model requires the following inputs: GHI, solar zenith angle, atmospheric pressure, and dew (wet-bulb) temperature.<jupyter_code># Estimation of the DIRINT model
dirint = pvlib.irradiance.dirint(
    ghi=df['ghi'],
    solar_zenith=solpos['apparent_zenith'],
    times=df.index,
    pressure=df['pressure']*100, # atmospheric pressure in Pascal 
    temp_dew=df['temp_dew']) # Dew temperature in Degree Celsius<jupyter_output><empty_output><jupyter_text>DIRINDEX modelThe DIRINDEX model (Pérez et al., 2002) was found to be one of the two best decomposition models in an extensive study by Gueymard and Ruiz-Arias (2016). The model is an enhancement of the DIRINT model by considering clearsky data.<jupyter_code># Retrieve clear-sky irradiance from CAMS McClear
clearsky, meta = pvlib.iotools.get_cams(
    latitude=location.latitude,
    longitude=location.longitude,
    altitude=location.altitude,
    start=df.index[0],
    end=df.index[-1],
    email=sodapro_username)

# Estimation of the DIRINDEX model
dirindex = pvlib.irradiance.dirindex(
    ghi=df['ghi'],
    ghi_clearsky=clearsky['ghi_clear'],
    dni_clearsky=clearsky['dni_clear'],
    zenith=solpos['zenith'],
    times=df.index,
    pressure=df['pressure']*100, # atmospheric pressure in Pascal 
    temp_dew=df['temp_dew']) # Dew temperature in Degree Celsius<jupyter_output><empty_output><jupyter_text>Erbs modelThe Erbs model (Erbs et al., 1982) is a simple and very popular model that only requires GHI and the solar zenith angle. This model is also available in pvlib and the implementation is equally straight-forward:<jupyter_code>erbs = pvlib.irradiance.erbs(
    ghi=df['ghi'],
    zenith=solpos['apparent_zenith'],
    datetime_or_doy=df.index)

erbs.columns<jupyter_output><empty_output><jupyter_text>Comparison and evaluation of decomposition modelsOnce the 4 models are estimated, we can evaluate the results of the models against the actual DNI measurements available. For example, below we visualize the results of the models a single day and a week in July:<jupyter_code># Plot of models for a single day
date = '2021-07-18'

disc.loc[date, 'dni'].plot(label='DISC', alpha=0.7)
dirint.loc[date].plot(label='DIRINT', alpha=0.7)
dirindex.loc[date].plot(label='DIRINDEX', alpha=0.7)
erbs.loc[date, 'dni'].plot(label='Erbs', alpha=0.7)
df.loc[date, 'dni'].plot(c='k', label='Reference')

plt.ylim(0, 1200)
plt.ylabel('Direct Normal Irradiance [W/m$^2$]')
plt.xlabel('Time [HH:MM]')
plt.title('Comparison of decomposition models\n(July 18th 2021)')
plt.legend(loc='upper left')
# Plot of models for several days on an hourly basis
disc.loc['20210709':'20210718', 'dni'].resample('1h').mean().plot(label='DISC', alpha=0.5) # format of date YYYYMMDD
dirint.loc['20210709':'20210718'].resample('1h').mean().plot(label='DIRINT', alpha=0.5)
dirindex.loc['20210709':'20210718'].resample('1h').mean().plot(label='DIRINDEX', alpha=0.5)
erbs.loc['20210709':'20210718', 'dni'].resample('1h').mean().plot(label='Erbs', alpha=0.5)
df.loc['20210709':'20210718', 'dni'].resample('1h').mean().plot(label='Reference')

plt.ylim(0,1200)
plt.ylabel('Direct Normal Irradiance [W/m$^2$]')
plt.title('Comparison of decomposition on an hourly basis\n (9th-18th July 2021)')
plt.legend(ncol=2, loc='upper center')
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>The correlation between the obtained models and the measured observations can be also presented in **scatter plots**:<jupyter_code># Create multiple plots
fig, axes = plt.subplots(ncols=4, figsize=(15, 4), facecolor='w', edgecolor='k', sharey=True)

# Plot
axes[0].scatter(df['dni'], disc['dni'], s=4, c='k')
axes[1].scatter(df['dni'], dirint, s=4, c='k')
axes[2].scatter(df['dni'], dirindex, s=4, c='k')
axes[3].scatter(df['dni'], erbs['dni'], s=4, c='k')

# Add characteristics to each subplot in a loop
axes[0].set_ylabel('Modelled DNI [W/m$^2$]')
axes[0].set_ylim(0, 1200)

for ax, model in zip(axes, ['DISC', 'DIRINT', 'DIRINDEX', 'ERBS']):
    ax.plot([0,1200], [0,1200], c='r')
    ax.set_xlim(0, 1200)
    ax.set_title(model)
    ax.set_xlabel('Observed DNI [W/m$^2$]')

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>The evaluation cannot be only based on a visual inspection. For solar resource assessment, there are a number of common metrics to evaluate the error and performance, which we will use below. These metrics are the mean bias error (MBE) and the root mean square error (RMSE). The MBE is the average error representing a systematic error to under- or overestimate and the RMSE is a measure of the dispersion of the deviations.$$ MBE = \frac{1}{n} \sum_{i=1}^{n} M_i - O_i $$$$ RMSE = \sqrt{\frac{1}{n}\sum_{i=1}^{n}{(M_i - O_i)^2}}$$where $n$ is the size of the sample or dataset, $M_i$ is each of the predicted values and $O_i$ are the observed or true values. In Python, each of these evaluation metrics can be defined as function.<jupyter_code># Mean Bias Error
def mbe(predictions, observed):
    return (predictions - observed).mean()

# Root Mean Square Error 
def rmse(predictions, observed):
    return np.sqrt(((predictions - observed)**2).mean())<jupyter_output><empty_output><jupyter_text>For the comparison of evaluation metrics, only measurement for periods where the solar elevation is greater than 5 degrees are used:<jupyter_code>df = df[solpos['elevation']>5]
erbs = erbs[solpos['elevation']>5]
disc = disc[solpos['elevation']>5]
dirint = dirint[solpos['elevation']>5]
dirindex = dirindex[solpos['elevation']>5]<jupyter_output><empty_output><jupyter_text>For the DISC method:<jupyter_code>mbe_disc = mbe(disc['dni'], df['dni'])
rmse_disc = rmse(disc['dni'], df['dni'])

print("DISC Decomposition Method:")
print("MBE " + "%.2f" % mbe_disc + " W/m\u00b2")
print("RMSE " + "%.2f" % rmse_disc + " W/m\u00b2")<jupyter_output>DISC Decomposition Method:
MBE 41.68 W/m²
RMSE 90.00 W/m²<jupyter_text>For the DIRINT method:<jupyter_code>mbe_dirint = mbe(dirint, df['dni'])
rmse_dirint = rmse(dirint, df['dni'])

print("DIRINT Decomposition Method:")
print("MBE " + "%.2f" % mbe_dirint + " W/m\u00b2")
print("RMSE " + "%.2f" % rmse_dirint + " W/m\u00b2")<jupyter_output>DIRINT Decomposition Method:
MBE 18.26 W/m²
RMSE 76.04 W/m²<jupyter_text>For the DIRINDEX method:<jupyter_code>mbe_dirindex = mbe(dirindex, df['dni'])
rmse_dirindex = rmse(dirindex, df['dni'])

print("DIRINDEX Decomposition Method:")
print("MBE " + "%.2f" % mbe_dirindex + " W/m\u00b2")
print("RMSE " + "%.2f" % rmse_dirindex + " W/m\u00b2")<jupyter_output>DIRINDEX Decomposition Method:
MBE -3.45 W/m²
RMSE 68.51 W/m²<jupyter_text>For the Erbs method:<jupyter_code>mbe_erbs = mbe(erbs['dni'], df['dni'])
rmse_erbs = rmse(erbs['dni'], df['dni'])

print("Erbs Decomposition Method:")
print("MBE " + "%.2f" % mbe_erbs + " W/m\u00b2")
print("RMSE " + "%.2f" % rmse_erbs + " W/m\u00b2")<jupyter_output>Erbs Decomposition Method:
MBE 14.06 W/m²
RMSE 81.79 W/m²
{'fraction_non_alphanumeric': 0.11453566621803499, 'fraction_numerical': 0.025975773889636608, 'mean_word_length': 5.892857142857143, 'pattern_counts': {'":': 0, '<': 82, '<?xml version=': 0, '>': 87, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '627186', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5099, 'n_tokens_neox': 4701, 'n_words': 1588, 'path': 'assessingsolar/notebooks/decomposition_models.ipynb', 'repo_name': 'AssessingSolar/Assessing-Solar', 'size': '14897'}