From 3054b6fce2e10c51ee861b620dda826351a5024f Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Mon, 2 Jun 2025 14:13:47 -0600 Subject: [PATCH 1/8] Create manual_corrections.csv --- metadata/manual_corrections.csv | 46 +++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 metadata/manual_corrections.csv diff --git a/metadata/manual_corrections.csv b/metadata/manual_corrections.csv new file mode 100644 index 0000000..bb74c83 --- /dev/null +++ b/metadata/manual_corrections.csv @@ -0,0 +1,46 @@ +start,end,columns,action,comment +"2015-01-01","2015-01-18","sun_diffuse;sun_total","remove" +"2015-01-01","2015-02-18","DNI;DHI","remove","30.1.15 er configfilen med de rigtige korrektionskonstanter til pyranometre sat i drift. Programmet læser ikke configfilen. Det problem er rettet af Lars Kokholm 18/2 2015" +"2015-01-01","2015-03-11 15:04","GHI","remove","11/3-15 kl 15:04:50: korrektionsfaktorerne for globalstrålingen f(CM11) og diffus stråling på vandret med skyggering (CM11) er ført ind i config-filen og cRIO er genstartet." +"2015-05-07 14:30","2015-05-13 19:20","GHI","remove","GHI was taken down for calibration (corresponds to Elsa's log file)" +"2017-09-06 08:45","2017-09-06 10:15","GHI","remove","GHI - Elsa's log states - 6/9 2017. CM11 til globalstråling er flyttet op på trackeren" +"2018-02-07 06:27","2018-02-07 16:00","DHI","remove","DHI pyranometer shows odd behavior (no diffuse on cloudy day which is clearly wrong)" +"2018-02-20 15:30","2018-02-20 16:30","DHI","remove","DHI went to zero" +"2018-03-08","2018-03-22","DHI","remove","DHI pyranometer stopped working (seems to be continuously decreasing starting March 8th 2018)" +"2018-03-21","2018-06-21 16:35","DHI","factor:8.34/8.59","DHI pyranometer 128764 (broken) was replaced by 128762 (sensitivity was first updated later)" +"2018-06-21 12:20","2018-06-21 15:10","DHI","remove","DHI pyranometer (128767) was recalibrated on June 21st 2018" +"2018-05-03","2018-06-08","DNI","remove","DNI was out for calibration in May and early June 2018" +"2018-07-20 18:30","2018-07-20 18:45","DNI","remove","Adjustment of the DNI sensor" +"2019-05-27 09","2025-03-31","sun_total;sun_diffuse","remove","SPN1 was taken down" +"2017-12-21 00","2018-01-11 09","sun_total;sun_diffuse","remove","Stuck at high value" +"2018-05-03","2018-05-13","weather","remove","WHICH ONE IS CORRECT" +"2018-05-03","2018-05-23","weather","remove","WHICH ONE IS CORRECT" +"2018-11-06 09:30","2018-11-06 12:00","weather","remove","Weather station down for check of rain sensor - rain sensor 'enabled' - Elsa" +"2018-05-23","2018-05-23 07:30","weather","remove","Weather station returned from calibration and had weird values" +"2018-08-13","2018-11-06","precipitation","remove","Rain sensor had been disabled in its configuration" +"2021-01-04 11:48","2021-01-27 16:15","DNI","remove","Instrument calibration" +"2021-01-04 12:46","2021-01-28 11:34","DHI","remove","Instrument calibration" +"2021-01-28 11:25","2021-02-22 17:25","GHI","remove","Instrument calibration" +"2023-01-25 14:56","2023-02-03 11:13","LWD","remove","Pyrgeometer was taken down for maitenance, specifically replacement of the bubble level. Dates correspond to log book entries" +"2023-03-28 04:00","2023-03-28 08:38","GHI","remove","Strange behavior of GHI pyranometer (measured much too low on clearsky day)" +"2024-03-07 13:47","2024-03-07 14:40","DNI","remove","DNI issue - check log book for other instruments?" +"2024-04-24 08:53","2024-04-24 14:24","GHI;DHI;DNI;LWD","remove","Change tracker location" +"2024-05-11 00:00","2024-05-16 23:59","GHI","remove","Humidity issues in GHI sensors" +"2024-08-17 05:00","2024-08-17 06:40","GHI;DHI;DNI","remove","Unknown deviation between GHI/GHI_calc - probably one sensor has dew" +"2024-09-19 05:00","2024-09-19 07:00","GHI;DHI;DNI","remove","DNI measurement too low - probably dew - clear-sky day" +"2024-09-21 05:00","2024-09-21 07:05","GHI;DHI;DNI","remove","DNI measurement too low - probably dew - clear-sky day" +"2015-01-01","2017-09-20","LWD","remove","Pyrgeometer measuring LWD was first correctly installed in September 2017" +"2017-11-06 18:00","2017-11-07 06:00","LWD","remove","sharp drop to negative valeus" +"2017-11-13 18:00","2017-11-14 06:00","LWD","remove","sharp drop to negative valeus" +"2017-10-31 05:36","2017-10-31 05:36","LWD","remove","sharp drop to negative valeus" +"2018-03-21 13:15","2018-03-21 13:32","LWD","remove","sharp increase" +"2018-11-07 00:00","2018-11-15 13:32","LWD","remove","sharp increase" +"2021-01-04 00:00","2021-02-23 00:00","LWD","remove","Out of service due to calibration" +"2025-03-27 12:00","2025-03-31 23:59","LWD","remove","Sensor was disconnected" +"2024-01-31 13:36","2024-01-31 13:37","LWD","remove","sharp increase" +"2024-03-07 14:37","2024-03-07 14:40","LWD","remove","sharp increase" +"2024-06-17 11:10","2024-06-17 11:26","LWD","remove","sharp increase" +"2024-03-29 17:43","2024-05-22 13:00","LWD","remove","Disconnected sensor cable to the LWD sensor" +"2025-03-27 13:15","2025-03-31 23:59","GHI;DHI;DNI;LWD","remove","Disconnected sensors, start of new time series." + +"2021-10-04 18:00","2021-10-16 00:00","DNI","remove","DNI is zero when it should not be" \ No newline at end of file From 98c26c432638436b5ec88048e5c9151ec3a56c18 Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Mon, 2 Jun 2025 14:13:51 -0600 Subject: [PATCH 2/8] Create parse_dtu_climate_station_data.py --- scripts/parse_dtu_climate_station_data.py | 373 ++++++++++++++++++++++ 1 file changed, 373 insertions(+) create mode 100644 scripts/parse_dtu_climate_station_data.py diff --git a/scripts/parse_dtu_climate_station_data.py b/scripts/parse_dtu_climate_station_data.py new file mode 100644 index 0000000..6b5513e --- /dev/null +++ b/scripts/parse_dtu_climate_station_data.py @@ -0,0 +1,373 @@ + +import glob +import pandas as pd +import pvlib +import pvanalytics +import os +import numpy as np +import warnings +import matplotlib.pyplot as plt + +data_path = 'C:/Users/arajen/OneDrive - Danmarks Tekniske Universitet/Skrivebord/Climate station data/' + +# %% To do +# When were the partly shaded pyranometers taken down? Make a list +# What is LWD really? Is it the upwelling or? + +# %% + +parameter_dict = { + 'Mod2Ch1': 'LWD', + 'Mod2Ch2': 'DHI_shadow_band', + 'Mod2Ch3': 'DHI', + 'Mod2Ch5': 'DNI', + 'Mod2Ch6': 'GHI', + 'Mod2Ch7': 'GTI_270_90', + 'Mod2Ch8': 'GTI_0_90', + 'Mod2Ch9': 'GTI_90_90', + 'Mod2Ch10': 'GTI_180_90', + 'Mod2Ch11': 'GTI_180_45', + 'RTD0_Pyrheliometer': 'DNI_temperature', + 'RTD2': 'LWD_temperature', +} + +weather_parameters = [ + 'wind_dir_min', 'wind_dir_avg', 'wind_dir_max', + 'wind_speed_min', 'wind_speed_avg', 'wind_speed_max', + 'air_temperature', 'relative_humidity', 'air_pressure', + 'rain_accumulation', 'rain_duration', 'rain_intensity', + 'hail_accumulation', 'hail_duration', 'hail_intensity', +] + + +precipitation_parameters = [ + 'rain_duration', 'rain_intensity', 'rain_accumulation', + 'hail_accumulation', 'hail_duration', 'hail_intensity', +] + +channels_not_in_use = [ + 'Mod2Ch4', # Not in use according to report + 'Mod3Ch1', 'Mod3Ch2', 'Mod3Ch3', 'Mod3Ch4', 'Mod3Ch5', 'Mod3Ch6', 'Mod3Ch7', 'Mod3Ch8', + 'Mod3Ch9', 'Mod3Ch10', 'Mod3Ch11', 'Mod3Ch12', 'Mod3Ch13', 'Mod3Ch14', 'Mod3Ch15', 'Mod3Ch16', + 'Mod2Ch12', 'Mod2Ch13', 'Mod2Ch14', 'Mod2Ch15', 'Mod2Ch16', +] + +half_dome_sensors = [f"Mod1Ch{i}" for i in range(1, 17)] + +irradiance_cols = ['GHI', 'DHI', 'DNI', 'sun_diffuse', 'sun_total'] + +gti_columns = ['GTI_270_90', 'GTI_0_90', 'GTI_90_90', 'GTI_180_90', 'GTI_180_45'] + +# %% Timestamp parsing + + +def date_parser(date_str): + """Parse timestamps in mixed formats.""" + if isinstance(date_str, (int, float)): + return pd.NaT + if date_str[2] == '-': # e.g., 04-01-2014 23:59 + format = '%m-%d-%Y %H:%M' + elif date_str[2] == ':': # e.g., 00:00:01 08/26/2014 + format = '%H:%M:%S %m/%d/%Y' + else: + raise UserWarning(f"Timestamp {date_str} does not match any format") + return pd.to_datetime(date_str, format=format) + + +# %% Sort filenames according to date + +def get_file_date(filename): + return pd.to_datetime(filename[-12:-4], format='%m-%d-%y') + + +filenames = glob.glob(data_path + '*') +filenames = [f for f in filenames if get_file_date(f).year != 2014] +file_dates = sorted([get_file_date(f) for f in filenames]) +filenames = sorted(filenames, key=lambda s: file_dates.index(get_file_date(s))) + +complete_dates = pd.date_range(min(file_dates), max(file_dates), freq='1d') +missing_dates = sorted(set(complete_dates).difference(file_dates)) + +print(pd.Series(missing_dates).dt.strftime('%Y-%m').value_counts()) + +# %% Read in data +dfs = [] +for f in filenames: + basename = os.path.basename(f) + file_date = get_file_date(f) + with warnings.catch_warnings(record=True) as w: + # Cause all warnings to always be triggered. + warnings.simplefilter("always") + + dfi = pd.read_csv(f, sep='\t', skip_blank_lines=True, + on_bad_lines='warn', na_values=[-9999]) + + dfi['file_date'] = file_date + + if len(w) == 1: # pd.read_csv engine choice specific + warning_messages = str(w[0].message).strip().split('\n') + warning_lines = [int(warning_messages[0].split()[2].strip(':')) + for m in warning_messages] + print(f"{basename.rstrip('.txt')} Warning lines: {len(warning_lines)}") + if len(w) > 1: + raise ValueError("Multiple warning types") + + dfs.append(dfi) + +# %% +df = pd.concat(dfs, axis=0) +df['Time(utc)'] = df['Time(utc)'].apply(lambda s: date_parser(s)) +df = df.set_index('Time(utc)') +df = df[df.index.notna()] +df = df.rename(columns=parameter_dict) +df.index = df.index.round('min') +# drop 2 instances of duplicated indexes (due to rounding) +df = df[~df.index.duplicated()] +# remove few incorrect time stamps (7 instances found in 2015 files) +df = df[df.index != '1904-01-01'] +# Ensure consistent frequency (add missing timestamps) +df = df.asfreq('1min') + + +# %% Drop unused/unnecessary columns +df = df.drop(columns=channels_not_in_use) +df = df.drop(columns=half_dome_sensors) +df = df.drop(columns=gti_columns) + +columns_drop = [ + 'wind_dir_min.1', # only available for part of January 2015 + 'Mod1ch16', # incorrect naming for January 2015 + # SOLYS 2 / sun sensor parameters + 'Q1_Sun_intensity', 'Q2_Sun_intensity', 'Q3_Sun_intensity', 'Q4_Sun_intensity', + 'Solys_instrument_status', 'Total_Sun_intensity', 'Function', + # Temperature parameters + 'RTD1', # seems to be just noise + 'RTD3', # value always stuck at 310.8 + 'RTD3_Cabinet', # values seem to be incorrect +] + +df = df.drop(columns=columns_drop) + +# %% Alternative NAN values + +nan_1290_parameters = ['sun_total', 'sun_diffuse', 'DNI_temperature', 'LWD_temperature'] + +for p in nan_1290_parameters: + df[p] = df[p].replace(-1290, np.nan) + +# %% Calculate solar position +location = pvlib.location.Location(latitude=55.7906, longitude=12.5251) +solpos = location.get_solarposition(df.index) + +df['GHI_calc'] = ( + df['DHI'].clip(lower=0) + + (df['DNI'] * np.cos(np.deg2rad(solpos['apparent_zenith']))).clip(lower=0)) + +df['dni_extra'] = pvlib.irradiance.get_extra_radiation(df.index) +df['ghi_extra'] = df['dni_extra'] * (np.cos(np.deg2rad(solpos['apparent_zenith']))).clip(lower=0) + +# %% + +df_qc = df.copy() + +# %% Apply manual corrections + +# Handle offset +# Calculate missing component + +qc_path = 'C:/github/dtu_solar_station/metadata/' +qc_filename = 'manual.csv' + + +def read_qc_file(filename): + """Read manual corrections file.""" + qci = pd.read_csv(filename, comment='#') + qci['start'] = pd.to_datetime(qci['start'], format='mixed') + qci['end'] = pd.to_datetime(qci['end'], format='mixed') + return qci + + +filename = qc_path + qc_filename +qc = read_qc_file(filename) + +for _, row in qc.iterrows(): + columns = row['columns'].split(';') + if ',' in row['columns']: + raise ValueError(f"Incorrectly formatted 'columns': {row['columns']}") + if row['start'] > row['end']: + raise ValueError(f"Start ({row['start']}) cannot be after end ({row['end']})") + # Apply same action to the pyrheliometer body temperature as the output signal + columns = columns + ['DNI_temperature'] if 'DNI' in columns else columns + columns = columns + ['LWD_temperature'] if 'LWD' in columns else columns + + if row['action'] == 'remove': + df_qc.loc[row['start']: row['end'], columns] = np.nan + + +# %% Plot all columns + +for c in df_qc.columns: + if c == 'Solys_gps_time': + continue + plt.figure() + df_qc.loc[df.index.year == 2015, c].plot() + plt.title(c) + plt.show() + +# %% +year = 2015 + +parameters = ['GHI', 'DHI', 'DNI'] + +df_qc.loc[df_qc.index.year==year, parameters].plot(ylim=(-10,1200), subplots=True, sharex=True) +#df_qc.loc['2018-01-01':, parameters].plot(ylim=(-200,600), subplots=True, sharex=True) + + +# %% BSRN checks of main measurements + +erl_dhi = pvanalytics.quality.irradiance.check_dhi_limits_qcrad( + dhi=df_qc['sun_diffuse'], solar_zenith=solpos['apparent_zenith'], + dni_extra=df_qc['dni_extra'], limits='extreme') +ppl_dhi = pvanalytics.quality.irradiance.check_dhi_limits_qcrad( + dhi=df_qc['sun_diffuse'], solar_zenith=solpos['apparent_zenith'], + dni_extra=df_qc['dni_extra'], limits='physical') + +erl_ghi = pvanalytics.quality.irradiance.check_ghi_limits_qcrad( + ghi=df_qc['sun_total'], solar_zenith=solpos['apparent_zenith'], + dni_extra=df_qc['dni_extra'], limits='extreme') + +ppl_spn1_ghi = pvanalytics.quality.irradiance.check_ghi_limits_qcrad( + ghi=df_qc['sun_total'], solar_zenith=solpos['apparent_zenith'], + dni_extra=df_qc['dni_extra'], limits='physical') + + +# %% +ghi_ppl, dhi_ppl, dni_ppl = pvanalytics.quality.irradiance.check_irradiance_limits_qcrad( + solar_zenith=solpos['apparent_zenith'], + dni_extra=df_qc['dni_extra'], + ghi=df_qc['GHI'], + dhi=df_qc['DHI'], + dni=df_qc['DNI'], + limits='extreme', +) + +# %% + +df_sub = df_qc['2021-09-25 00:00': '2021-10-20 00:00'] +df_sub = df_qc['2021-09-25 00:00': '2021-10-06 00:00'] +df_sub = df_qc['2021-10-04 00:00': '2021-10-06 00:00'] +axes = df_sub[['GHI','DHI','DNI']].plot(subplots=True, sharex=True) + + + +plot_measured_vs_calculated_ghi(df_sub['GHI'], df_sub['DHI'], df_sub['DNI'], + solpos.loc[df_sub.index, 'apparent_zenith'], s=1, alpha=0.5) +# _ = plot_diffuse_fraction_vs_clearness_index( +# ghi=df_sub['GHI'], dhi=df_sub['DHI'], solar_zenith=solpos.loc[df_sub.index, 'apparent_zenith'], +# dni_extra=df_sub['dni_extra'], mask=df_sub['GHI']>50, s=0.5, alpha=0.2, xlim=(None, None), ylim=(None, None)) + +# %% + + +# %% + +def plot_diffuse_fraction_vs_clearness_index(ghi, dhi, solar_zenith, dni_extra, mask=None, + c=None, xlim=(0, 2), ylim=(0, 1.5), **kwargs): + mask = np.ones(len(ghi)).astype(bool) if mask is None else mask + c = c if c is None else c[mask] + + fig, ax = plt.subplots() + plt.scatter( + x = ghi[mask] / (dni_extra[mask] * np.cos(np.deg2rad(solar_zenith[mask])).clip(lower=0)), + y = dhi[mask] / ghi[mask], + c=c, **kwargs) + + ax.set_xlabel('Diffuse fraction (DHI/GHI) [-]') + ax.set_ylabel('Clearness index (GHI/ETH) [-]') + ax.set_xlim(xlim) + ax.set_ylim(ylim) + plt.show() + return fig, ax + + +def plot_measured_vs_calculated_ghi(ghi, dhi, dni, solar_zenith, mask=None, c=None, + xlim=(-10, 1200), ylim=(-10, 1200), **kwargs): + mask = np.ones(len(ghi)).astype(bool) if mask is None else mask + c = c if c is None else c[mask] + + fig, ax = plt.subplots() + plt.scatter( + x = ghi[mask], + y = dhi[mask].clip(lower=0) + \ + (dni[mask] * np.cos(np.deg2rad(solar_zenith[mask])).clip(lower=0)).clip(lower=0), + c=c, **kwargs) + + ax.set_xlabel('Measured GHI [W/m$^2$]') + ax.set_ylabel('Calculated GHI [W/m$^2$]') + ax.plot([0, 1200], [0, 1200], c='r', alpha=0.8) + ax.set_xlim(xlim) + ax.set_ylim(ylim) + plt.show() + return fig, ax + +_ = plot_diffuse_fraction_vs_clearness_index( + ghi=df_qc['GHI'], dhi=df_qc['DHI'], solar_zenith=solpos['apparent_zenith'], + dni_extra=df_qc['dni_extra'], mask=df_qc['GHI']>50, s=0.5, alpha=0.2) + +# %% SPN1 check + +dfp = df_qc[df_qc['sun_total'] > 50] + +fig, ax = plt.subplots() +plt.scatter( + x=dfp['sun_total'] / dfp['ghi_extra'], + y=dfp['sun_diffuse'] / dfp['sun_total'], + c=solpos.loc[dfp.index, 'apparent_elevation'], + s=0.5, alpha=0.2) + +ax.set_xlabel('sun_total / ghi_extra') +ax.set_ylabel('sun_diffuse / sun_total') +# ax.set_xlim(0, 2) +# ax.set_ylim(-0.1, 1.2) + +# %% SPN1 BSRN checks +consistent_components_spn1, diffuse_ratio_limit_spn1 = \ + pvanalytics.quality.irradiance.check_irradiance_consistency_qcrad( + solar_zenith=solpos['apparent_zenith'], + ghi=df_qc['sun_total'], + dhi=df_qc['sun_diffuse'], + dni=np.nan, + outside_domain=True) + +erl_spn1_dhi = pvanalytics.quality.irradiance.check_dhi_limits_qcrad( + dhi=df_qc['sun_diffuse'], solar_zenith=solpos['apparent_zenith'], + dni_extra=df_qc['dni_extra'], limits='extreme') +ppl_spn1_dhi = pvanalytics.quality.irradiance.check_dhi_limits_qcrad( + dhi=df_qc['sun_diffuse'], solar_zenith=solpos['apparent_zenith'], + dni_extra=df_qc['dni_extra'], limits='physical') + +erl_spn1_ghi = pvanalytics.quality.irradiance.check_ghi_limits_qcrad( + ghi=df_qc['sun_total'], solar_zenith=solpos['apparent_zenith'], + dni_extra=df_qc['dni_extra'], limits='extreme') + +ppl_spn1_ghi = pvanalytics.quality.irradiance.check_ghi_limits_qcrad( + ghi=df_qc['sun_total'], solar_zenith=solpos['apparent_zenith'], + dni_extra=df_qc['dni_extra'], limits='physical') + +# %% +df_qc.loc['2024-06-17 11:09':'2024-06-17 11:27', 'LWD'].plot(ylim=(-200,600)) + + +# %% + + +df_qc.loc[(df.index.year==2015), irradiance_cols].plot( + subplots=True, sharex=True, figsize=(8,8)) + +# %% +# 2018, 2021 was an ugly year + +df_qc[df_qc.index.year==202].plot.scatter( + x='GHI', y='sun_total', xlim=[-10, 1200], ylim=[-10, 1200], + s=0.1, alpha=0.4) +plt.plot([0, 1200], [0, 1200], c='r', alpha=0.5) From 5f365d91cb7a28f0c4458b654c41c0076ec44c21 Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Tue, 3 Jun 2025 10:30:25 -0600 Subject: [PATCH 3/8] Update parse_dtu_climate_station_data.py --- scripts/parse_dtu_climate_station_data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/parse_dtu_climate_station_data.py b/scripts/parse_dtu_climate_station_data.py index 6b5513e..9bbbb7f 100644 --- a/scripts/parse_dtu_climate_station_data.py +++ b/scripts/parse_dtu_climate_station_data.py @@ -176,7 +176,7 @@ def get_file_date(filename): # Calculate missing component qc_path = 'C:/github/dtu_solar_station/metadata/' -qc_filename = 'manual.csv' +qc_filename = 'manual_corrections.csv' def read_qc_file(filename): From 121f488a698af62a42805710bf9604e3c11fe4aa Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Wed, 4 Jun 2025 19:08:41 -0600 Subject: [PATCH 4/8] Update manual_corrections.csv Update based on investigations done using multiplots and Elsa's old logbook --- metadata/manual_corrections.csv | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/metadata/manual_corrections.csv b/metadata/manual_corrections.csv index bb74c83..954160a 100644 --- a/metadata/manual_corrections.csv +++ b/metadata/manual_corrections.csv @@ -6,15 +6,16 @@ start,end,columns,action,comment "2017-09-06 08:45","2017-09-06 10:15","GHI","remove","GHI - Elsa's log states - 6/9 2017. CM11 til globalstråling er flyttet op på trackeren" "2018-02-07 06:27","2018-02-07 16:00","DHI","remove","DHI pyranometer shows odd behavior (no diffuse on cloudy day which is clearly wrong)" "2018-02-20 15:30","2018-02-20 16:30","DHI","remove","DHI went to zero" -"2018-03-08","2018-03-22","DHI","remove","DHI pyranometer stopped working (seems to be continuously decreasing starting March 8th 2018)" +"2018-03-08 00:00","2018-03-21 23:59","DHI","remove","DHI pyranometer stopped working (seems to be continuously decreasing starting March 8th 2018)" "2018-03-21","2018-06-21 16:35","DHI","factor:8.34/8.59","DHI pyranometer 128764 (broken) was replaced by 128762 (sensitivity was first updated later)" +"'2015-03-11 00:00","2015-05-13 17:20","GHI","factor:4.98/5.065","Then the new calibration coefficient (4.98) for the GHI pyranometer was put in place on May 13 2015, teh closure equation was much better. Therefore, this coefficient is backfilled for the first five months of the year." "2018-06-21 12:20","2018-06-21 15:10","DHI","remove","DHI pyranometer (128767) was recalibrated on June 21st 2018" -"2018-05-03","2018-06-08","DNI","remove","DNI was out for calibration in May and early June 2018" +"2018-05-03 00:00","2018-06-07 11:00","DNI","remove","DNI was out for calibration in May and early June 2018" "2018-07-20 18:30","2018-07-20 18:45","DNI","remove","Adjustment of the DNI sensor" -"2019-05-27 09","2025-03-31","sun_total;sun_diffuse","remove","SPN1 was taken down" -"2017-12-21 00","2018-01-11 09","sun_total;sun_diffuse","remove","Stuck at high value" -"2018-05-03","2018-05-13","weather","remove","WHICH ONE IS CORRECT" -"2018-05-03","2018-05-23","weather","remove","WHICH ONE IS CORRECT" +"2019-03-25 00:00","2019-05-27 09:00","sun_total;sun_diffuse","remove","Logbook notes that there was dew inside the SPN1." +"2019-05-27 09:00","2025-03-31","sun_total;sun_diffuse","remove","SPN1 was taken down" +"2017-12-21 00:00","2018-01-11 09","sun_total;sun_diffuse","remove","Stuck at high value" +"2018-05-03 00:00","2018-05-23","weather","remove","The weather station was down for calibration." "2018-11-06 09:30","2018-11-06 12:00","weather","remove","Weather station down for check of rain sensor - rain sensor 'enabled' - Elsa" "2018-05-23","2018-05-23 07:30","weather","remove","Weather station returned from calibration and had weird values" "2018-08-13","2018-11-06","precipitation","remove","Rain sensor had been disabled in its configuration" @@ -29,7 +30,7 @@ start,end,columns,action,comment "2024-08-17 05:00","2024-08-17 06:40","GHI;DHI;DNI","remove","Unknown deviation between GHI/GHI_calc - probably one sensor has dew" "2024-09-19 05:00","2024-09-19 07:00","GHI;DHI;DNI","remove","DNI measurement too low - probably dew - clear-sky day" "2024-09-21 05:00","2024-09-21 07:05","GHI;DHI;DNI","remove","DNI measurement too low - probably dew - clear-sky day" -"2015-01-01","2017-09-20","LWD","remove","Pyrgeometer measuring LWD was first correctly installed in September 2017" +"2015-01-01 00:00","2017-09-20 23:59","LWD","remove","Pyrgeometer measuring LWD was first correctly installed in September 2017" "2017-11-06 18:00","2017-11-07 06:00","LWD","remove","sharp drop to negative valeus" "2017-11-13 18:00","2017-11-14 06:00","LWD","remove","sharp drop to negative valeus" "2017-10-31 05:36","2017-10-31 05:36","LWD","remove","sharp drop to negative valeus" @@ -42,5 +43,11 @@ start,end,columns,action,comment "2024-06-17 11:10","2024-06-17 11:26","LWD","remove","sharp increase" "2024-03-29 17:43","2024-05-22 13:00","LWD","remove","Disconnected sensor cable to the LWD sensor" "2025-03-27 13:15","2025-03-31 23:59","GHI;DHI;DNI;LWD","remove","Disconnected sensors, start of new time series." - -"2021-10-04 18:00","2021-10-16 00:00","DNI","remove","DNI is zero when it should not be" \ No newline at end of file +"2021-10-04 18:00","2021-10-16 00:00","DNI","remove","DNI is zero when it should not be" +"2016-04-14 06:15","2016-04-14 07:26","DHI","remove","Unresonable low DHI value, could be due or other obstacle" +"2018-01-07 07:40","2018-01-07 12:00","DHI","remove","Believed to be icing or snow cover of the diffuse sensor" +"2018-01-08 07:40","2018-01-08 12:00","DHI","remove","Believed to be icing or snow cover of the diffuse sensor" +"2017-09-14 11:24","2017-09-14 12:22","DHI","remove","DHI sensor output dropped to zero for about an hour." +"2017-09-15 07:25","2017-09-15 12:40","DHI","remove","DHI sensor output dropped to zero for about an hour." +"2017-09-16 00:00","2017-09-16 13:15","DHI","remove","DHI sensor output is too low" +"2024-02-01 07:49","2024-03-07 15:30","GHI","replace","Ventilator heater was on but the ventilator itself was not powered" \ No newline at end of file From 844a161021a72d65efdc4e19f08b8c3acdd8e452 Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Wed, 4 Jun 2025 20:09:54 -0600 Subject: [PATCH 5/8] Add solar angle cut offs --- metadata/manual_corrections.csv | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/metadata/manual_corrections.csv b/metadata/manual_corrections.csv index 954160a..7dd6148 100644 --- a/metadata/manual_corrections.csv +++ b/metadata/manual_corrections.csv @@ -50,4 +50,5 @@ start,end,columns,action,comment "2017-09-14 11:24","2017-09-14 12:22","DHI","remove","DHI sensor output dropped to zero for about an hour." "2017-09-15 07:25","2017-09-15 12:40","DHI","remove","DHI sensor output dropped to zero for about an hour." "2017-09-16 00:00","2017-09-16 13:15","DHI","remove","DHI sensor output is too low" -"2024-02-01 07:49","2024-03-07 15:30","GHI","replace","Ventilator heater was on but the ventilator itself was not powered" \ No newline at end of file +"2024-02-01 07:49","2024-03-07 15:30","GHI","replace","Ventilator heater was on but the ventilator itself was not powered" +"2024-01-19 00:00","2025-03-31 23:59","Azimuth;Azimuth_adjust;Zenith;Zenith_adjust","remove","Angles returned from the Solys2 were zero or missing for this period." From ba1dd442618dedfe43e8b1e908a21b707a4fc464 Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Wed, 4 Jun 2025 21:44:24 -0600 Subject: [PATCH 6/8] Additional manual corrections from going into details in plots --- metadata/manual_corrections.csv | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/metadata/manual_corrections.csv b/metadata/manual_corrections.csv index 7dd6148..c6449e4 100644 --- a/metadata/manual_corrections.csv +++ b/metadata/manual_corrections.csv @@ -52,3 +52,11 @@ start,end,columns,action,comment "2017-09-16 00:00","2017-09-16 13:15","DHI","remove","DHI sensor output is too low" "2024-02-01 07:49","2024-03-07 15:30","GHI","replace","Ventilator heater was on but the ventilator itself was not powered" "2024-01-19 00:00","2025-03-31 23:59","Azimuth;Azimuth_adjust;Zenith;Zenith_adjust","remove","Angles returned from the Solys2 were zero or missing for this period." +"2017-02-24 00:00","2017-02-24 12:00","GHI","remove","Likely snow or icing of the pyranometer. GHI has strange signal and is lower than calculated GHI." +"2018-03-20 00:00","2018-03-20 10:00","GHI","remove","Kn is higher than Kt. Perfect clear sky day, so issue must be with GHI which also has a strange profile." +"2018-03-21 00:00","2018-03-21 14:30","DHI","remove","DHI was near zero, unrealistic. Perhas a loose wire." +"2020-11-27 07:30","2020-11-27 10:00","DNI","remove","GHI_calc is lower than GHI and DNI seems to deviate from clear sky profile." +"2020-11-28 08:00","2020-11-28 11:00","GHI","remove","GHI does not follow clear sky profile whereas DHI and DNI does." +"2021-12-24 12:00","2021-12-25 23:59","GHI","remove","GHI is too low for a clear-sky day. Perhaps covered in layer of snow or ice." +"2022-12-16 00:00","2022-12-16 23:59","GHI","remove","GHI is too low for a clear-sky day. Perhaps covered in layer of snow or ice." +"2023-03-02 06:00","2023-03-02 09:00","GHI","remove","GHI has dew strikes." From e8357412967478208c6d7806d91fe17ee4bf75f5 Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Wed, 4 Jun 2025 22:14:03 -0600 Subject: [PATCH 7/8] Update parse_dtu_climate_station_data.py --- scripts/parse_dtu_climate_station_data.py | 182 ++-------------------- 1 file changed, 17 insertions(+), 165 deletions(-) diff --git a/scripts/parse_dtu_climate_station_data.py b/scripts/parse_dtu_climate_station_data.py index 9bbbb7f..0d5871a 100644 --- a/scripts/parse_dtu_climate_station_data.py +++ b/scripts/parse_dtu_climate_station_data.py @@ -88,7 +88,7 @@ def get_file_date(filename): complete_dates = pd.date_range(min(file_dates), max(file_dates), freq='1d') missing_dates = sorted(set(complete_dates).difference(file_dates)) -print(pd.Series(missing_dates).dt.strftime('%Y-%m').value_counts()) +print("Missing days\n", pd.Series(missing_dates).dt.strftime('%Y-%m').value_counts()) # %% Read in data dfs = [] @@ -129,6 +129,11 @@ def get_file_date(filename): df = df.asfreq('1min') +# Note from Elsa's logbook: +# 11-5-2018: P1-P16 samt kupler er pillet ned. Dog sidder pyranometeret Mod1Ch7 stadig og måler +# med en helt lukket kuppel. +df.loc['2018-05-11':, half_dome_sensors] = np.nan + # %% Drop unused/unnecessary columns df = df.drop(columns=channels_not_in_use) df = df.drop(columns=half_dome_sensors) @@ -166,7 +171,7 @@ def get_file_date(filename): df['dni_extra'] = pvlib.irradiance.get_extra_radiation(df.index) df['ghi_extra'] = df['dni_extra'] * (np.cos(np.deg2rad(solpos['apparent_zenith']))).clip(lower=0) -# %% +# %% Copy dataframe df_qc = df.copy() @@ -204,170 +209,17 @@ def read_qc_file(filename): df_qc.loc[row['start']: row['end'], columns] = np.nan -# %% Plot all columns - -for c in df_qc.columns: - if c == 'Solys_gps_time': - continue - plt.figure() - df_qc.loc[df.index.year == 2015, c].plot() - plt.title(c) - plt.show() - -# %% -year = 2015 - -parameters = ['GHI', 'DHI', 'DNI'] - -df_qc.loc[df_qc.index.year==year, parameters].plot(ylim=(-10,1200), subplots=True, sharex=True) -#df_qc.loc['2018-01-01':, parameters].plot(ylim=(-200,600), subplots=True, sharex=True) - - -# %% BSRN checks of main measurements - -erl_dhi = pvanalytics.quality.irradiance.check_dhi_limits_qcrad( - dhi=df_qc['sun_diffuse'], solar_zenith=solpos['apparent_zenith'], - dni_extra=df_qc['dni_extra'], limits='extreme') -ppl_dhi = pvanalytics.quality.irradiance.check_dhi_limits_qcrad( - dhi=df_qc['sun_diffuse'], solar_zenith=solpos['apparent_zenith'], - dni_extra=df_qc['dni_extra'], limits='physical') - -erl_ghi = pvanalytics.quality.irradiance.check_ghi_limits_qcrad( - ghi=df_qc['sun_total'], solar_zenith=solpos['apparent_zenith'], - dni_extra=df_qc['dni_extra'], limits='extreme') - -ppl_spn1_ghi = pvanalytics.quality.irradiance.check_ghi_limits_qcrad( - ghi=df_qc['sun_total'], solar_zenith=solpos['apparent_zenith'], - dni_extra=df_qc['dni_extra'], limits='physical') - - -# %% -ghi_ppl, dhi_ppl, dni_ppl = pvanalytics.quality.irradiance.check_irradiance_limits_qcrad( - solar_zenith=solpos['apparent_zenith'], - dni_extra=df_qc['dni_extra'], - ghi=df_qc['GHI'], - dhi=df_qc['DHI'], - dni=df_qc['DNI'], - limits='extreme', -) - -# %% +# These +df_qc.loc['2018-03-21 00:00': '2018-06-21 16:35', 'DHI'] = 8.34/8.59 * \ + df_qc.loc['2018-03-21 00:00': '2018-06-21 16:35', 'DHI'] -df_sub = df_qc['2021-09-25 00:00': '2021-10-20 00:00'] -df_sub = df_qc['2021-09-25 00:00': '2021-10-06 00:00'] -df_sub = df_qc['2021-10-04 00:00': '2021-10-06 00:00'] -axes = df_sub[['GHI','DHI','DNI']].plot(subplots=True, sharex=True) +df_qc.loc['2015-03-11 00:00': '2015-05-13 17:20', 'GHI'] = 4.98/5.065 * \ + df_qc.loc['2015-03-11 00:00': '2015-05-13 17:20', 'GHI'] +# Important to recalculate +df_qc['GHI_calc'] = ( + df_qc['DHI'].clip(lower=0) + + (df_qc['DNI'] * np.cos(np.deg2rad(solpos['apparent_zenith']))).clip(lower=0)) -plot_measured_vs_calculated_ghi(df_sub['GHI'], df_sub['DHI'], df_sub['DNI'], - solpos.loc[df_sub.index, 'apparent_zenith'], s=1, alpha=0.5) -# _ = plot_diffuse_fraction_vs_clearness_index( -# ghi=df_sub['GHI'], dhi=df_sub['DHI'], solar_zenith=solpos.loc[df_sub.index, 'apparent_zenith'], -# dni_extra=df_sub['dni_extra'], mask=df_sub['GHI']>50, s=0.5, alpha=0.2, xlim=(None, None), ylim=(None, None)) - -# %% - - -# %% - -def plot_diffuse_fraction_vs_clearness_index(ghi, dhi, solar_zenith, dni_extra, mask=None, - c=None, xlim=(0, 2), ylim=(0, 1.5), **kwargs): - mask = np.ones(len(ghi)).astype(bool) if mask is None else mask - c = c if c is None else c[mask] - - fig, ax = plt.subplots() - plt.scatter( - x = ghi[mask] / (dni_extra[mask] * np.cos(np.deg2rad(solar_zenith[mask])).clip(lower=0)), - y = dhi[mask] / ghi[mask], - c=c, **kwargs) - - ax.set_xlabel('Diffuse fraction (DHI/GHI) [-]') - ax.set_ylabel('Clearness index (GHI/ETH) [-]') - ax.set_xlim(xlim) - ax.set_ylim(ylim) - plt.show() - return fig, ax - - -def plot_measured_vs_calculated_ghi(ghi, dhi, dni, solar_zenith, mask=None, c=None, - xlim=(-10, 1200), ylim=(-10, 1200), **kwargs): - mask = np.ones(len(ghi)).astype(bool) if mask is None else mask - c = c if c is None else c[mask] - - fig, ax = plt.subplots() - plt.scatter( - x = ghi[mask], - y = dhi[mask].clip(lower=0) + \ - (dni[mask] * np.cos(np.deg2rad(solar_zenith[mask])).clip(lower=0)).clip(lower=0), - c=c, **kwargs) - - ax.set_xlabel('Measured GHI [W/m$^2$]') - ax.set_ylabel('Calculated GHI [W/m$^2$]') - ax.plot([0, 1200], [0, 1200], c='r', alpha=0.8) - ax.set_xlim(xlim) - ax.set_ylim(ylim) - plt.show() - return fig, ax - -_ = plot_diffuse_fraction_vs_clearness_index( - ghi=df_qc['GHI'], dhi=df_qc['DHI'], solar_zenith=solpos['apparent_zenith'], - dni_extra=df_qc['dni_extra'], mask=df_qc['GHI']>50, s=0.5, alpha=0.2) - -# %% SPN1 check - -dfp = df_qc[df_qc['sun_total'] > 50] - -fig, ax = plt.subplots() -plt.scatter( - x=dfp['sun_total'] / dfp['ghi_extra'], - y=dfp['sun_diffuse'] / dfp['sun_total'], - c=solpos.loc[dfp.index, 'apparent_elevation'], - s=0.5, alpha=0.2) - -ax.set_xlabel('sun_total / ghi_extra') -ax.set_ylabel('sun_diffuse / sun_total') -# ax.set_xlim(0, 2) -# ax.set_ylim(-0.1, 1.2) - -# %% SPN1 BSRN checks -consistent_components_spn1, diffuse_ratio_limit_spn1 = \ - pvanalytics.quality.irradiance.check_irradiance_consistency_qcrad( - solar_zenith=solpos['apparent_zenith'], - ghi=df_qc['sun_total'], - dhi=df_qc['sun_diffuse'], - dni=np.nan, - outside_domain=True) - -erl_spn1_dhi = pvanalytics.quality.irradiance.check_dhi_limits_qcrad( - dhi=df_qc['sun_diffuse'], solar_zenith=solpos['apparent_zenith'], - dni_extra=df_qc['dni_extra'], limits='extreme') -ppl_spn1_dhi = pvanalytics.quality.irradiance.check_dhi_limits_qcrad( - dhi=df_qc['sun_diffuse'], solar_zenith=solpos['apparent_zenith'], - dni_extra=df_qc['dni_extra'], limits='physical') - -erl_spn1_ghi = pvanalytics.quality.irradiance.check_ghi_limits_qcrad( - ghi=df_qc['sun_total'], solar_zenith=solpos['apparent_zenith'], - dni_extra=df_qc['dni_extra'], limits='extreme') - -ppl_spn1_ghi = pvanalytics.quality.irradiance.check_ghi_limits_qcrad( - ghi=df_qc['sun_total'], solar_zenith=solpos['apparent_zenith'], - dni_extra=df_qc['dni_extra'], limits='physical') - -# %% -df_qc.loc['2024-06-17 11:09':'2024-06-17 11:27', 'LWD'].plot(ylim=(-200,600)) - - -# %% - - -df_qc.loc[(df.index.year==2015), irradiance_cols].plot( - subplots=True, sharex=True, figsize=(8,8)) - -# %% -# 2018, 2021 was an ugly year - -df_qc[df_qc.index.year==202].plot.scatter( - x='GHI', y='sun_total', xlim=[-10, 1200], ylim=[-10, 1200], - s=0.1, alpha=0.4) -plt.plot([0, 1200], [0, 1200], c='r', alpha=0.5) +# df_qc.to_pickle('quality_controlled_dtu_data.pkl') From 63fbd9209341afd1258fd738b9bdd28f8265763a Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Fri, 21 Nov 2025 00:18:04 +0100 Subject: [PATCH 8/8] Update parse_dtu_climate_station_data.py --- scripts/parse_dtu_climate_station_data.py | 88 +++++++++++++++++++++++ 1 file changed, 88 insertions(+) diff --git a/scripts/parse_dtu_climate_station_data.py b/scripts/parse_dtu_climate_station_data.py index 0d5871a..f18542a 100644 --- a/scripts/parse_dtu_climate_station_data.py +++ b/scripts/parse_dtu_climate_station_data.py @@ -223,3 +223,91 @@ def read_qc_file(filename): + (df_qc['DNI'] * np.cos(np.deg2rad(solpos['apparent_zenith']))).clip(lower=0)) # df_qc.to_pickle('quality_controlled_dtu_data.pkl') + +# %% Export to Task 16 format + +meta = { + 'station_code': 'LYN', + 'latitude': '55.79065', + 'longitude': '12.52509', + 'altitude': '40', + 'timezone': '1', +} + +df_t16 = df_qc.copy() + +df_t16 = df_t16.rename(columns={'DHI': 'DIF'}) + +full_index = pd.date_range( + start=df_t16.index.min().round('min').replace(month=1, day=1), + end=df_t16.index.max().round('min').replace(month=12, day=31, hour=23, minute=59), + freq='1min') + +df_t16 = df_t16.reindex(full_index) + +# %% +# create individual columns containing the year, month,day, hour, and minute +# based on the dataframe index +df_t16['Year'] = df_t16.index.year +df_t16['Month'] = df_t16.index.month +df_t16['Day'] = df_t16.index.day +df_t16['Hour'] = df_t16.index.hour +df_t16['Minute'] = df_t16.index.minute + + +empty_cols = [ + 'GHIcalc', 'Elev', 'Azim', 'Kc', 'usable', 'flagPPLDIF', 'flagERLDIF', + 'flagPPLDNI', 'flagERLDNI', 'flagPPLGHI', 'flagERLGHI', 'flag3highSZA', + 'flag3lowSZA', 'flagKt', 'flagKKt', 'flagKn', 'flagKnKt', 'flagKhighSZA', + 'flagKlowSZA', 'flagTracker', 'flagManual'] + +df_t16[empty_cols] = np.nan + +columns_in_specific_order = [ + 'Year', 'Month', 'Day', 'Hour', 'Minute', + 'GHI', 'DNI', 'DIF', 'GHIcalc', + 'Elev', 'Azim', 'Kc', 'usable', + 'flagPPLDIF', 'flagERLDIF', 'flagPPLDNI', 'flagERLDNI', + 'flagPPLGHI', 'flagERLGHI', 'flag3highSZA', 'flag3lowSZA', + 'flagKt', 'flagKKt', 'flagKn', 'flagKnKt', + 'flagKhighSZA', 'flagKlowSZA', 'flagTracker', 'flagManual' +] + +df_t16 = df_t16[columns_in_specific_order] + +# %% Generate header lines + +header = [ + f"# stationcode {meta['station_code']},,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + f"# latitude deg N {meta['latitude']},,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + f"# longitude deg E {meta['longitude']},,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + f"# altitude in m amsl {meta['altitude']},,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + f"# timezone offset from UTC in hours {meta['timezone']},,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + "# time stamp (Year Month Day Hour Minute) reference is in UTC and refers to the end of the period,,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + "# missing data points in GHI DNI and DIF are noted with -999,,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + "# GHI is the global horizontal irradiance in W/m2,,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + "# DNI is the direct normal irradiance in W/m2,,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + "# DIF is the diffuse horizontal irradiance in W/m2,,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + "# GHIcalc is the calculated GHI from DNI and DIF,,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + "# Elev is the solar elevation in deg,,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + "# Azim is the solar azimuth angle in deg N,,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + "# Kc is the clearsky index calculated with CAMS mcclear with GHIcalc/GHI McClear,,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + "# usable is the validity of a data point with 1 being valid and 0 being not usable,,,,,,,,,,,,,,,,,,,,,,,,,,,,,", + "# flag values from various tests: 1 means the data failed the test. 0 means the test was passed. -999 means the test domain was not met,,,,,,,,,,,,,,,,,,,,,,,,,,,,,", +] + +header = '\n'.join(header) + '\n' + + +# %% Save the data in files based on the operation year + +output_path = 'C:/users/arajen/downloads/' + +for year, group in df_t16.groupby('Year'): + # Open the file for writing + with open(output_path + f"{meta['station_code']}_{year}.csv", "w", newline='') as f: + # Write the header lines at the top of the CSV file + f.write(header) + # Write the DataFrame content to the CSV file + group.to_csv(f, index=False, header=True) + print(f"Saved {meta['station_code']}_{year}.csv")