5. 2013 HOPE

processed with pyrnet-0.2.16

The PyrNet was setup for calibration during HOPE-Melpitz on the Melpitz measurement field. Cross-calibration is done versus reference observations from the TROPOS MObile RaDiation ObseRvatory (MORDOR) station.

5.1. Imports

#|dropcode
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 pvlib import clearsky
from pvlib.location import Location

import pyrnet.pyrnet

5.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. Here we generate the calibration.json file for the processing to l1b:

box_numbers = np.arange(1,101)
calibrations = {f"{bn:03d}":[7,7] for bn in box_numbers}
calibjson = {"2000-01-01": calibrations}
with open("pyrnet_calib_prep.json","w") as txt:
    json.dump(calibjson, txt)

Within pyrnet_config.json:

{"file_calibration" : "pyrnet_calib_prep.json"}

Workflow for preparation

  1. Prepare pyrnet_config_calibration_prep.json with contributors metadata and the dummy calibration config file.

  2. $ pyrnet process l1a -c pyrnet_config.json raw_data/*.bin l1a/

  3. $ pyrnet process l1b_network -c pyrnet_config.json l1a/*.nc l1b_network/

5.3. Configuration

Set local data paths and lookup metadata.

pf_mordor = "mordor/CR1000_Radiation_corrected.dat"
pf_pyrnet = "l1b_network/{date:%Y-%m-%d}_P1D_pyrnet_hope-melpitz_n000l1bf1s.c01.nc"

loc = Location(51.525175, 12.91648, altitude=90) # Melpitz

stations = np.arange(1,101)
# 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 )

5.3.1. Load reference MORDOR data

#|dropcode
#|dropout

df = pd.read_csv(
    pf_mordor,
    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')
cs_mask = clearsky.detect_clearsky(
    dst['TP2_Wm2'].values,
    cs['ghi'],
    times=pd.to_datetime(dst.time.values)
)    

dst = dst.assign({"cs_mask":("time", cs_mask)})

    
mordor = dst.copy()
mordor = mordor.drop_duplicates("time", keep="last")
mordor
<xarray.Dataset>
Dimensions:  (time: 11394)
Coordinates:

  • time (time) datetime64[ns] 2013-09-19T12:47:00 … 2013-09-27T10:40:00 Data variables: TP2_Wm2 (time) float64 597.0 713.0 762.6 764.3 … 401.6 482.4 577.5 646.0 cs_mask (time) bool False False False False … False False False False

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='.')
ax.plot(cs.index,cs["ghi"])
[<matplotlib.lines.Line2D at 0x7f0eee1cb970>]

png

5.3.2. Load PyrNet Data

#|dropcode
#|dropout
dates = pd.date_range(
    pd.to_datetime(mordor.time.values[0].astype("datetime64[D]")),
    pd.to_datetime(mordor.time.values[-1].astype("datetime64[D]")),
    freq="1d"
)

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:          (station: 50, maintenancetime: 4, time: 5299)
Coordinates:

  • station (station) int64 2 7 14 16 18 22 23 … 89 90 92 95 96 100

  • maintenancetime (maintenancetime) datetime64[ns] 2013-09-19T11:55:00 ……

  • time (time) datetime64[ns] 2013-09-19T05:59:00 … 2013-09-27… Data variables: ghi (time, station) float64 0.0001123 nan … 0.0005218 szen (time, station) float64 79.97 79.97 79.97 … 79.96 79.96 Attributes: (12/31) Conventions: CF-1.10, ACDD-1.3 title: TROPOS pyranometer network (PyrNet) observatio… history: 2024-11-13T10:53:01: 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: 2013-09-19T00:00:00 time_coverage_end: 2013-09-19T23:59:59 time_coverage_duration: P0DT23H59M59S time_coverage_resolution: P0DT0H0M1S site: [‘’, ‘’, ‘’, ‘’, ‘’, ‘’, ‘’, ‘’, ‘’, ‘’, ‘’, ‘…

5.4. Calibration

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

TODO: Revise versus 2023 EU version.

Cloudy sky treatment is applied.

5.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 $

# # Set flux values to nan if no pyranometer is installed.
# pyr.ghi.values = pyr.ghi.where(mainmask).values

# # convert to measured voltage
# pyr.ghi.values = pyr.ghi.values * 7 * 1e-6

# # Step 1, select data
# pyr = pyr.where(pyr.szen<80, drop=True)
# pyr.ghi.values = pyr.ghi.where(pyr.ghi>0.033/300.).values

5.4.2. Step 2

Interpolate reference to PyrNet samples and combine to a single Dataset

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

# Calibration datasets for main and extra pyranometer
Cds_main = xr.Dataset(
    data_vars={
        'reference_Wm2': ('time', mordor['TP2_Wm2'].data),
        'pyrnet_V': (('time','station'), pyr['ghi'].data)
    },
    coords= {
        "time": pyr.time,
        "station": pyr.station
    }
)

5.4.3. Step 3

Remove outliers from series using xarray grouping and apply function. The following functions removes outliers (deviation more than 2% according to ISO 9847) from a selected group. This step involves calculating calibration series and the integration of one hour intervals to smooth out high variable situation, which would break the calibration even when time synchronization is slightly off. Also this gets rid of some random shading events like birds / chimney / rods in line of sigth, which would affect calibration otherwise. We following ISO 9847 5.4.1.1 equation (2) here.

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_V'] / x['reference_Wm2']
    # integrated series 
    ix = x.integrate('time')
    M = ix['pyrnet_V'] / ix['reference_Wm2']
    
    while np.any(np.abs(C-M) > 0.02*M):
        #calculate as long there are outliers deviating more than 2 percent
        x = x.where(np.abs(C-M) < 0.02*M)
        C = x['pyrnet_V'] / x['reference_Wm2']
        #integrated series 
        ix = x.integrate('time')
        M = ix['pyrnet_V'] / 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)


5.4.4. Step 4

The series of measured voltage and irradiance is now without outliers. So we use equation 1 again to calculate from this reduced series the calibration factor for the instant samples.

C_main = 1e6*Cds_main['pyrnet_V'] / Cds_main['reference_Wm2']
C_main.values[C_main.values<6]=np.nan
C_main.values[C_main.values>8]=np.nan

5.4.5. Step 5

We just found the Calibration factor to be the mean of the reduced calibration factor series and the uncertainty to be the standard deviation of this reduced series. Steo 3, 4 and 5 are done for every pyranometer seperate.

C_main_mean = C_main.mean(dim='time',skipna=True)
C_main_std = C_main.std(dim='time',skipna=True)

5.5. Results

#|dropcode
fig, ax = plt.subplots(1,1, figsize=(10,5))
ax.set_title("Main Pyranometer")
ax.plot(C_main.time, C_main, ls ="", marker='.')
ax.set_xlabel("Date")
ax.set_ylabel("Calibration factor (uV / Wm-2)")
ax.grid(True)
fig.show()

plt.figure()
fig, ax = plt.subplots(1,1, figsize=(10,5))
ax.set_title("Main Pyranometer")
ax.plot(pyr.szen.interp_like(C_main), C_main, ls ="", marker='.')
ax.set_xlabel("solar zenith angle (deg)")
ax.set_ylabel("Calibration factor (uV / Wm-2)")
ax.grid(True)
fig.show()



png

<Figure size 640x480 with 0 Axes>

png

calibration_new = {}
print(f"Box:    Main       ,     Extra  ")
for box in C_main_mean.station:
    Cm = float(C_main_mean.sel(station=box).values)
    Um = float(C_main_std.sel(station=box).values)
    
    calibration_new.update({
        f"{box:03d}": [np.round(Cm,2), None]
    })
    print(f"{box:3d}: {Cm:.2f} +- {Um:.3f} , {None}")

calibjson = {"2013-09-01": calibration_new}
with open("pyrnet_calib_new.json","w") as txt:
    json.dump(calibjson, txt)
Box:    Main       ,     Extra  
  2: 7.48 +- 0.300 , None
  7: 7.39 +- 0.297 , None
 14: 7.43 +- 0.366 , None
 16: 7.60 +- 0.278 , None
 18: 7.15 +- 0.460 , None
 22: 7.42 +- 0.206 , None
 23: 7.48 +- 0.240 , None
 28: 7.51 +- 0.188 , None
 29: 7.09 +- 0.297 , None
 32: 7.52 +- 0.441 , None
 35: 7.18 +- 0.374 , None
 37: 7.49 +- 0.342 , None
 38: 7.07 +- 0.456 , None
 40: 7.43 +- 0.233 , None
 42: 7.45 +- 0.240 , None
 43: 7.25 +- 0.316 , None
 48: 7.14 +- 0.432 , None
 49: 7.17 +- 0.381 , None
 51: 7.22 +- 0.283 , None
 53: 7.50 +- 0.268 , None
 54: 6.95 +- 0.503 , None
 56: 6.58 +- 0.289 , None
 58: 7.47 +- 0.266 , None
 60: 7.01 +- 0.486 , None
 63: 7.27 +- 0.375 , None
 65: 7.12 +- 0.397 , None
 66: 6.84 +- 0.386 , None
 67: 6.87 +- 0.510 , None
 68: 6.78 +- 0.330 , None
 69: 7.16 +- 0.330 , None
 70: 6.79 +- 0.399 , None
 71: 7.35 +- 0.344 , None
 73: 7.08 +- 0.495 , None
 74: 7.20 +- 0.383 , None
 75: 6.62 +- 0.270 , None
 77: 7.34 +- 0.236 , None
 78: 6.62 +- 0.255 , None
 79: 7.38 +- 0.291 , None
 80: 6.86 +- 0.303 , None
 81: 7.27 +- 0.253 , None
 85: 7.21 +- 0.391 , None
 86: 7.13 +- 0.267 , None
 87: 7.22 +- 0.079 , None
 88: 7.05 +- 0.266 , None
 89: 7.23 +- 0.413 , None
 90: 7.28 +- 0.367 , None
 92: 6.95 +- 0.330 , None
 95: 7.60 +- 0.111 , None
 96: 7.29 +- 0.399 , None
100: 7.00 +- 0.267 , None