4. 2014 IOPRAO¶

processed with pyrnet-0.2.16

The PyrNet was setup for calibration in a dense array on the Melpitz measurement field from 2015-05-06 to 2015-05-11. Cross-calibration is done versus reference observations from the TROPOS MObile RaDiation ObseRvatory (MORDOR) station.

4.1. Imports¶

#|dropcode
from IPython.display import display, Latex
import os
import xarray as xr
import pandas as pd
import numpy as np
import datetime as dt
import matplotlib.pyplot as plt
import jstyleson as json

from scipy.optimize import differential_evolution

from pvlib import clearsky
from pvlib.location import Location

import pyrnet.pyrnet
import pyrnet.utils

4.2. Prepare PyrNet data¶

For calibration preparation the PyrNet data is processed to level l1b using a calibration factor of 7 (uV W-1 m2) for all pyranometers with the pyrnet process l1b tool. This is done to unify the conversion to sensor voltage during calibration and not run into valid_range limits for netcdf encoding.

Before running this notebook, new absolute calibration factors have to be determined with calibration_melcol.ipynb

4.3. Configuration¶

Set local data paths and lookup metadata.

pf_mordor = "mordor/{date:%Y/%m/%Y-%m-%d}_Radiation.dat"
pf_bsrn = "bsrn/{date:%Y-%m}_bsrn.tab"
pf_pyrnet = "l1b_network/{date:%Y-%m-%d}_P1D_pyrnet_ioprao_n000l1bf1s.c01.nc"
dates = pd.date_range("2014-06-04","2014-07-18")
# period with lots of clear sky situations
dates = pd.date_range("2014-06-06","2014-06-09")
stations = np.arange(1,101)

loc = Location(52.21, 14.122, altitude=125) # Lindenberg

# lookup which box contains actually a pyranometer/ extra pyranometer
mainmask = [] 
for box in stations:
    _, serials, _, _ = pyrnet.pyrnet.meta_lookup(dates[0],box=box)
    mainmask.append( True if len(serials[0])>0 else False )

4.4. Calibration¶

The calibration follows the ISO 9847:1992 - Solar energy — Calibration of field pyranometers by comparison to a reference pyranometer. For selecting and masking the data.

TODO: Revise versus 2023 EU version.

Cloudy sky treatment is applied.

4.4.1. Step 1¶

Drop Night measures and low signal measures from pyranometer data. Since calibration without incoming radiation doesnt work.

This data is kept for calibration:

solar zenith angle < 80° ( as recommended in ISO 9847)
Measured Voltage > 0.033 V, e.g. ADC count is 0 or 1 of 1023 (drop the lowest ~1%)

Voltage measured ($V_m$) at the logger is the amplified Senor voltage ($V_S$) by a gain of 300.

$ V_m = 300 V_S$

As the uncalibrated flux measurements ($F_U$) are calibrated with a fixed factor of 7 uV W-1 m2:

$ V_s = 71e-6 F_U $

4.4.2. Step 2¶

Apply new absolute calibration from calibration_ioprao.ipynb

calib_new = pyrnet.utils.read_json("pyrnet_calib_new_bsrn.json")
calib_new = calib_new[list(calib_new.keys())[0]]
for station in pyr.station:
    pyr.ghi.sel(station=station).values /= calib_new[f"{station:03d}"][0] * 1e-6 

4.4.3. Step 3¶

Interpolate reference to PyrNet samples and combine to a single Dataset

# interpolate reference to PyrNet
mordor = mordor.interp(time=pyr.time).interpolate_na()
bsrn = bsrn.interp(time=pyr.time).interpolate_na()

# Calibration datasets for main and extra pyranometer
Cds_main = xr.Dataset(
    data_vars={
        'reference2_Wm2': ('time', mordor['TP2_Wm2'].data[bsrn.cs_mask]),
        'reference_Wm2': ('time', bsrn['SWD'].data[bsrn.cs_mask]),
        'pyrnet_Wm2': (('time','station'), pyr['ghi'].data[bsrn.cs_mask,:]),
        'szen': ('time',pyr.szen.mean("station",skipna=True).data[bsrn.cs_mask])
    },
    coords= {
        "time": pyr.time[bsrn.cs_mask],
        "station": pyr.station
    }
)

4.4.4. Step 4¶

Remove outliers from series using xarray grouping and apply function.

def remove_outliers(x):
    """
    x is an xarray dataset containing these variables:
    coords: 'time' - datetime64
    'pyrnet_V' - array - voltage measures of pyranometer
    'reference_Wm2' - array - measured irradiance of reference
    """

    # calculate calibration series for single samples
    C = x['pyrnet_Wm2'] / x['reference_Wm2']
    # integrated series 
    ix = x.integrate('time')
    M = ix['pyrnet_Wm2'] / ix['reference_Wm2']
    
    while np.any(np.abs(C-M) > 0.01*M):
        #calculate as long there are outliers deviating more than 2 percent
        x = x.where(np.abs(C-M) < 0.01*M)
        C = x['pyrnet_Wm2'] / x['reference_Wm2']
        #integrated series 
        ix = x.integrate('time')
        M = ix['pyrnet_Wm2'] / ix['reference_Wm2']
        
    #return the reduced dataset x
    return x

# remove outliers
Cds_main = Cds_main.groupby('time.hour').apply(remove_outliers)

# hourly mean
Cds_main = Cds_main.resample(time="1h").mean(skipna=True)

4.4.5. Step 5¶

The series of measured reference and pyrnet irradiance is now without outliers. Now we can fit a cubic function - (reference/pyrnet) vs. szen - to determine the cosine correction function for pyrnet.

fig,(ax,ax2) = plt.subplots(1,2)
p = ax.plot(Cds_main.szen,Cds_main.reference_Wm2/Cds_main.pyrnet_Wm2,ls='',marker='.')
p = ax2.plot(Cds_main.szen,Cds_main.reference2_Wm2/Cds_main.pyrnet_Wm2,ls='',marker='.')

ax.grid(True)
ax2.grid(True)

png

def cubic_function(x, a, b, c, d):
    return a * x ** 3 + b * x ** 2 + c * x + d

def objective_function(coefficients, x, y):
    a, b, c, d = coefficients
    y_pred = cubic_function(x, a, b, c, d)
    error = np.sum((y - y_pred) ** 2)
    return error
bounds = [(-4, -2), (3, 7), (-5, 0), (0, 3)]


ratio = Cds_main.reference_Wm2/Cds_main.pyrnet_Wm2
ratio = ratio.where(ratio<1.3)
ratio2 = Cds_main.reference2_Wm2/Cds_main.pyrnet_Wm2
ratio2 = ratio2.where(ratio2<1.3)

result = differential_evolution(
    objective_function,
    args=(np.cos(np.deg2rad(Cds_main.szen)), ratio),
    bounds=bounds,
    seed=1
)

result2 = differential_evolution(
    objective_function,
    args=(np.cos(np.deg2rad(Cds_main.szen)), ratio2),
    bounds=bounds,
    seed=1
)

print(result)
print(result2)

 message: Optimization terminated successfully.
 success: True
     fun: 0.6873391346566766
       x: [-3.327e+00  5.846e+00 -3.141e+00  1.489e+00]
     nit: 47
    nfev: 3010
     jac: [-4.976e-05  8.521e-05  1.756e-04 -3.683e-04]
 message: Optimization terminated successfully.
 success: True
     fun: 0.7129148950067137
       x: [-3.110e+00  5.538e+00 -3.000e+00  1.463e+00]
     nit: 51
    nfev: 3175
     jac: [-1.436e-04  7.412e-05 -1.138e-04  1.075e-04]

4.5. Results¶

print("Best cubic fit (BSRN):")
a3, a2, a1, a0 = result.x
display(Latex(
    rf"""
    {a3:+.3f}x^3{a2:+.3f}x^2{a1:+.3f}x{a0:+.3f}
    """
))
print("Best cubic fit (MORDOR):")
a32, a22, a12, a02 = result2.x
display(Latex(
    rf"""
    {a32:+.3f}x^3{a22:+.3f}x^2{a12:+.3f}x{a02:+.3f}
    """
))

Best cubic fit (BSRN):

-3.327x^3+5.846x^2-3.141x+1.489

Best cubic fit (MORDOR):

-3.110x^3+5.538x^2-3.000x+1.463

# Coefficients from other calibrations:

display(Latex(
    rf"""
    IOPRAO BSRN: {a3:+.3f}x^3{a2:+.3f}x^2{a1:+.3f}x{a0:+.3f}
    """
))
display(Latex(
    rf"""
    IOPRAO MORDOR: {a32:+.3f}x^3{a22:+.3f}x^2{a12:+.3f}x{a02:+.3f}
    """
))

b3, b2, b1, b0 = [-2.33, 4.4,-2.44, 1.38]
display(Latex(
    rf"""
    MelCol: {b3:+.3f}x^3{b2:+.3f}x^2{b1:+.3f}x{b0:+.3f}
    """
))


c3, c2, c1, c0 = [ -3.01, 5.59, -3.04, 1.45 ]
display(Latex(
    rf"""
    MetPVNet: {c3:+.3f}x^3{c2:+.3f}x^2{c1:+.3f}x{c0:+.3f}
    """
))
d3, d2, d1, d0 = [ -2.227, 4.366, -2.524, 1.385 ]
display(Latex(
    rf"""
    S2VSR: {d3:+.3f}x^3{d2:+.3f}x^2{d1:+.3f}x{d0:+.3f}
    """
))

IOPRAO BSRN: -3.327x^3+5.846x^2-3.141x+1.489

IOPRAO MORDOR: -3.110x^3+5.538x^2-3.000x+1.463

MelCol: -2.330x^3+4.400x^2-2.440x+1.380

MetPVNet: -3.010x^3+5.590x^2-3.040x+1.450

S2VSR: -2.227x^3+4.366x^2-2.524x+1.385

szen = np.arange(1,80)
mu0 = np.cos(np.deg2rad(szen))

fig,ax = plt.subplots(1,1)
p = ax.plot(Cds_main.szen,ratio,ls='',marker='.')
ax.plot(szen, a3*mu0**3 + a2*mu0**2 + a1*mu0 + a0,
         color='C0', label=f'best-cubic-fit (BSRN): {a3:+.3f}x^3{a2:+.3f}x^2{a1:+.3f}x{a0:+.3f}')
ax.plot(szen, a32*mu0**3 + a22*mu0**2 + a12*mu0 + a02,
         color='C1', label=f'best-cubic-fit (MORDOR): {a32:+.3f}x^3{a22:+.3f}x^2{a12:+.3f}x{a02:+.3f}')
ax.plot(szen, b3*mu0**3 + b2*mu0**2 + b1*mu0 + b0,
         color='C2',ls=':', label=f'MelCol: {b3:+.3f}x^3{b2:+.3f}x^2{b1:+.3f}x{b0:+.3f}')
ax.plot(szen, c3*mu0**3 + c2*mu0**2 + c1*mu0 + c0,
         color='C3',ls=':', label=f'MetPvNet: {c3:+.3f}x^3{c2:+.3f}x^2{c1:+.3f}x{c0:+.3f}')
ax.plot(szen, d3*mu0**3 + d2*mu0**2 + d1*mu0 + d0,
         color='C4',ls=':', label=f'S2VSR: {d3:+.3f}x^3{d2:+.3f}x^2{d1:+.3f}x{d0:+.3f}')

ax.legend()
ax.set_xlabel('Zenith angle [deg] ')
ax.set_ylabel('reference_ghi/pyr_ghi // correction_factor [-]')
ax.grid(True)

png

# dump to json
calibjson = {"2014-06-01": {"CC":list(result.x[::-1])}}
with open("pyrnet_calib_cc_new_bsrn.json","w") as txt:
    json.dump(calibjson, txt)

calibjson = {"2014-06-01": {"CC":list(result2.x[::-1])}}
with open("pyrnet_calib_cc_new_mordor.json","w") as txt:
    json.dump(calibjson, txt)

4. 2014 IOPRAO¶

4.1. Imports¶

4.2. Prepare PyrNet data¶

4.3. Configuration¶

4.3.1. Load reference MORDOR data¶

4.3.2. Load BSRN data¶

4.3.3. Load PyrNet Data¶