## 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. ### Imports ```python #|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 ``` ### 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 ### Configuration Set local data paths and lookup metadata. ```python 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 ) ``` #### Load reference MORDOR data The reference data of MORDOR is loaded, and clearsky is detected on daily basis using solis_simple clearsky model and the pvlib.clearsky.detect_clearsky function. ```python #|dropcode #|dropout new = True for i,date in enumerate(dates): fname = pf_mordor.format(date=date) if not os.path.exists(fname): continue df = pd.read_csv( fname, header=0, skiprows=[0,2,3], date_format="ISO8601", na_values=["NAN"], parse_dates=[0], index_col=0 ) dst = df.to_xarray().rename({"TIMESTAMP":"time"}) # drop not needed variables keep_vars = ['TP2_Wm2'] # global shortwave irradiance drop_vars = [v for v in dst if v not in keep_vars] dst = dst.drop_vars(drop_vars) dst = dst.resample(time="1min").mean(skipna=True) cs = loc.get_clearsky(pd.to_datetime(dst.time.values),model='simplified_solis') try: cs_mask = clearsky.detect_clearsky( dst['TP2_Wm2'].values, cs['ghi'], times=pd.to_datetime(dst.time.values) ) except: cs_mask = np.zeros(dst.time.size).astype(bool) dst = dst.assign({"cs_mask":("time", cs_mask)}) # merge if new: ds = dst.copy() new = False else: ds = xr.concat((ds,dst),dim='time', data_vars='minimal', coords='minimal', compat='override') mordor = ds.copy() mordor = mordor.drop_duplicates("time", keep="last") mordor ```

<xarray.Dataset>
Dimensions:  (time: 5719)
Coordinates:
  * time     (time) datetime64[ns] 2014-06-06 ... 2014-06-09T23:18:00
Data variables:
    TP2_Wm2  (time) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... nan nan nan nan nan 0.0
    cs_mask  (time) bool False False False False ... False False False False

```python fig,ax = plt.subplots(1,1) ax.plot(mordor.time,mordor.TP2_Wm2) ax.plot(mordor.time[mordor.cs_mask],mordor.TP2_Wm2[mordor.cs_mask],ls='',marker='.') ``` [] ![png](calibration_ioprao_cc_output_8_1.png) #### Load BSRN data ```python umonth = np.unique(dates.values.astype("datetime64[M]")) new = True for month in umonth: fname = pf_bsrn.format(date=pd.to_datetime(month)) if not os.path.exists(fname): continue df = pd.read_csv( fname, sep='\s+', header=None, skiprows=34, date_format="ISO8601", parse_dates=[0], index_col=0, names=["time","SWD"], usecols=[0,2] ) dst = df.to_xarray() dst.SWD.values[dst.SWD.values<0] = 0 dst.SWD.values = dst.SWD.values.astype(float) cs = loc.get_clearsky(pd.to_datetime(dst.time.values),model='simplified_solis') try: cs_mask = clearsky.detect_clearsky( dst['SWD'].values, cs['ghi'].values, times=pd.to_datetime(dst.time.values), window_length=60 ) except: cs_mask = np.zeros(dst.time.size).astype(bool) dst = dst.assign({"cs_mask":("time", cs_mask)}) # merge if new: ds = dst.copy() new = False else: ds = xr.concat((ds,dst),dim='time', data_vars='minimal', coords='minimal', compat='override') bsrn = ds.copy() bsrn = bsrn.drop_duplicates("time", keep="last") bsrn ```

<xarray.Dataset>
Dimensions:  (time: 43200)
Coordinates:
  * time     (time) datetime64[ns] 2014-06-01 ... 2014-06-30T23:59:00
Data variables:
    SWD      (time) float64 1.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0 0.0
    cs_mask  (time) bool False False False False ... False False False False

#### Load PyrNet Data ```python #|dropcode #|dropout for i,date in enumerate(dates): # read from thredds server dst = xr.open_dataset(pf_pyrnet.format(date=date)) # drop not needed variables keep_vars = ['ghi','szen'] drop_vars = [v for v in dst if v not in keep_vars] dst = dst.drop_vars(drop_vars) # unify time and station dimension to speed up merging date = dst.time.values[0].astype("datetime64[D]") timeidx = pd.date_range(date, date + np.timedelta64(1, 'D'), freq='1s', inclusive='left') dst = dst.interp(time=timeidx) dst = dst.reindex({"station": stations}) dst.ghi.values = dst.ghi.values * 7 * 1e-6 dst = dst.where(dst.szen<80, drop=True) dst.ghi.values = dst.ghi.where(dst.ghi>0.033/300.).values dst = dst.resample(time="1min").mean(skipna=True) # merge if i == 0: ds = dst.copy() else: ds = xr.concat((ds,dst),dim='time', data_vars='minimal', coords='minimal', compat='override') pyr = ds.copy() pyr ```

<xarray.Dataset>
Dimensions:          (time: 3319, station: 43, maintenancetime: 1)
Coordinates:
  * station          (station) int64 1 5 6 8 15 16 20 ... 90 91 92 93 95 98 99
  * maintenancetime  (maintenancetime) datetime64[ns] 2014-06-11T23:59:59
  * time             (time) datetime64[ns] 2014-06-06T04:08:00 ... 2014-06-09...
Data variables:
    ghi              (time, station) float64 0.0003242 0.0003215 ... nan nan
    szen             (time, station) float64 79.99 79.99 79.99 ... 79.99 nan nan
Attributes: (12/31)
    Conventions:               CF-1.10, ACDD-1.3
    title:                     TROPOS pyranometer network (PyrNet) observatio...
    history:                   2024-11-12T10:03:39: Merged level l1b by pyrne...
    institution:               Leibniz Institute for Tropospheric Research (T...
    source:                    TROPOS pyranometer network (PyrNet)
    references:                https://doi.org/10.5194/amt-9-1153-2016
    ...                        ...
    geospatial_lon_units:      degE
    time_coverage_start:       2014-06-06T00:00:00
    time_coverage_end:         2014-06-06T23:59:59
    time_coverage_duration:    P0DT23H59M59S
    time_coverage_resolution:  P0DT0H0M1S
    site:                      ['', '', '', '', '', '', '', '', '', '', '', '...

### Calibration The calibration follows the [ISO 9847:1992 - Solar energy — Calibration of field pyranometers by comparison to a reference pyranometer](https://archive.org/details/gov.in.is.iso.9847.1992). For selecting and masking the data. > TODO: Revise versus 2023 EU version. Cloudy sky treatment is applied. #### 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 = 7*1e-6* F_U $ #### Step 2 Apply new absolute calibration from `calibration_ioprao.ipynb` ```python 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 ``` #### Step 3 Interpolate reference to PyrNet samples and combine to a single Dataset ```python # 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 } ) ``` #### Step 4 Remove outliers from series using xarray grouping and apply function. ```python 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) ``` #### 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. ```python 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](calibration_ioprao_cc_output_22_0.png) ```python 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] ### Results ```python 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 ```python # 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 ```python 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](calibration_ioprao_cc_output_27_0.png) ```python # 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) ```