Skip to content

Commit

Permalink
Update SWE_postproc_RTKLIB_Neumayer.py
Browse files Browse the repository at this point in the history
  • Loading branch information
@lasteine
lasteine committed Jun 22, 2023
1 parent 9aaa0a8 commit d4f67ac
Showing 1 changed file with 2 additions and 163 deletions.
165 changes: 2 additions & 163 deletions SWE_postproc_RTKLIB_Neumayer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,11 +37,7 @@
rover = 'ReachM2_sladina-raw_' # 'NMER' or '3393' (old Emlid: 'ReachM2_sladina-raw_')
rover_name = 'NMER_original' # 'NMER' or 'NMER_original' or 'NMLR'
receiver = 'NMER' # 'NMER' or 'NMLB' or 'NMLR'
<<<<<<< HEAD
base_name = 'NMLB' # prefix of base rinex observation files, e.g. station name
=======
base_name = 'NMLB'
>>>>>>> 05d3630 (Optimized for reading only new data textfiles and attach them to the binary pickle.)
base = '3387' # rinex file name prefix for base receiver
nav = '3387' # navigation file name prefix for broadcast ephemerides files
sp3 = 'COD' # navigation file name prefix for precise ephemerides files
Expand All @@ -50,26 +46,16 @@
options_Leica = 'rtkpost_options_Ladina_Leica_statisch_multisystemfrequency_neumayer_900_15' # name of RTKLIB configuration file (.conf) for high-end receiver
options_Emlid = 'rtkpost_options_Ladina_Emlid_statisch_multisystemfrequency_neumayer_900_15' # name of RTKLIB configuration file (.conf) for low-cost receiver
ending = '' # file name suffix if needed: e.g., a variant of the processing '_eleambmask15', '_noglonass'
<<<<<<< HEAD
acc_y_lim = (-100, 2000) # y-axis limit for accumulation plots
delta_acc_y_lim = (-200, 1000) # y-axis limit for delta accumulation plots
swe_y_lim = (-100, 700) # y-axis limit for water equivalent plots
delta_swe_y_lim = (-100, 500) # y-axis limit for delta water equivalent plots
xlim_dates = dt.date(2021, 11, 26), dt.date(2023, 4, 1) # time series date limits to plot on x-axis
cal_date = '2022-07-24' # calibration date for snow density estimation
yy = '21' # initial year of observations
=======
acc_y_lim = (-200, 1600) # y-axis limit for accumulation plots
delta_acc_y_lim = (-400, 1000) # y-axis limit for delta accumulation plots
swe_y_lim = (-100, 700) # y-axis limit for water equivalent plots
delta_swe_y_lim = (-200, 600) # y-axis limit for delta water equivalent plots
xlim_dates = dt.date(2021, 11, 26), dt.date(2023, 1, 17) # time series date limits to plot on x-axis
offset_le = 28 # offset between leica and emlid rover antennas in mm w.e
yy = '21' # initial year of observations
>>>>>>> 05d3630 (Optimized for reading only new data textfiles and attach them to the binary pickle.)
save_plots = True # show (False) or save (True) plots
total_backup = False # copy (True) all new data to server for backup, else (False) do not copy
solplot_backup = True # copy (True) all new solution files and plots to server for backup, else (False) do not copy
solplot_backup = False # copy (True) all new solution files and plots to server for backup, else (False) do not copy


""" 0. Preprocess data """
Expand All @@ -87,11 +73,7 @@
# check available solution data to only further process new data (first mjd in 2021 here: 59544.0)
yy, start_mjd, start_mjd_emlid = f.get_sol_yeardoy(dest_path, resolution)

<<<<<<< HEAD
# # for proceeding manually without preprocessing: calculate start/end mjd using a given start/end date
=======
# # for proceedina manually without preprocessing: calculate start/end mjd using a given start/end date
>>>>>>> 05d3630 (Optimized for reading only new data textfiles and attach them to the binary pickle.)
# start_mjd, end_mjd = f.get_mjd_int(2021, 11, 12, 2023, 2, 28)
# start_mjd_emlid, end_mjd_emlid = start_mjd, end_mjd

Expand All @@ -109,26 +91,11 @@


''' 3. Filter and clean ENU solution data '''
<<<<<<< HEAD
fil_df_emlid, u_emlid, u_clean_emlid, swe_gnss_emlid, std_gnss_emlid, swe_gnss_daily_emlid, std_gnss_daily_emlid = f.filter_rtklib_solutions(
dest_path, 'NMER', resolution, df_enu=df_enu_emlid, ambiguity=1, threshold=3, window='D', ending=ending)

fil_df_leica, u_leica, u_clean_leica, swe_gnss_leica, std_gnss_leica, swe_gnss_daily_leica, std_gnss_daily_leica = f.filter_rtklib_solutions(
dest_path, 'NMLR', resolution, df_enu=df_enu_leica, ambiguity=1, threshold=3, window='D', ending=ending)
=======
# filter and clean ENU solution data (outlier filtering, median filtering, adjustments for observation mast heightening) and store results in pickle and .csv
fil_df_emlid, fil_emlid, fil_clean_emlid, m_emlid, s_emlid, jump_emlid, swe_gnss_emlid, swe_gnss_daily_emlid, std_gnss_daily_emlid = f.filter_rtklib_solutions(
dest_path, 'NMER', resolution, df_enu=df_enu_emlid, ambiguity=1, ti_set_swe2zero=12, threshold=3, window='D',
resample=False, resample_resolution='30min', ending=ending)

fil_df_leica, fil_leica, fil_clean_leica, m_leica, s_leica, jump_leica, swe_gnss_leica, swe_gnss_daily_leica, std_gnss_daily_leica = f.filter_rtklib_solutions(
dest_path, 'NMLR', resolution, df_enu=df_enu_leica, ambiguity=1, ti_set_swe2zero=12, threshold=3, window='D',
resample=False, resample_resolution='30min', ending=ending)
>>>>>>> 05d3630 (Optimized for reading only new data textfiles and attach them to the binary pickle.)

# correct offset leica to emlid
swe_gnss_emlid = swe_gnss_emlid - offset_le
swe_gnss_daily_emlid = swe_gnss_daily_emlid - offset_le


""" 4. Read reference sensors data """
Expand All @@ -145,12 +112,7 @@
leica_daily, emlid_daily, buoy_daily, poles_daily, laser_daily = f.resample_allobs(gnss_leica, gnss_emlid, buoy, poles, laser_filtered, interval='D')


<<<<<<< HEAD
""" 6. Calculate SWE differences, linear regressions & RMSE between GNSS refractometr and reference data """
=======
""" 6. Calculate differences, linear regressions, RMSE & MRB between GNSS and reference data """
# TODO: calculate correct values (with correct input data)
>>>>>>> 05d3630 (Optimized for reading only new data textfiles and attach them to the binary pickle.)
# calculate differences between reference data and GNSS (Leica/Emlid)
diffs_sh_daily, diffs_swe_daily = f.calculate_differences2gnss(emlid_daily, leica_daily, manual, laser_daily, buoy_daily, poles_daily)
diffs_sh_15min, diffs_swe_15min, laser_15min = f.calculate_differences2gnss_15min(gnss_emlid, gnss_leica, laser_filtered)
Expand All @@ -165,11 +127,7 @@
f.calculate_rmse(diffs_swe_daily, diffs_swe_15min, manual, laser_filtered, gnssir_acc=None, gnssir_acc_daily=None, res='SWE')


<<<<<<< HEAD
''' 7. Plot GNSS Refractometry results (SWE, ΔSWE, scatter) '''
=======
''' 7. Plot results (SWE, ΔSWE, scatter) '''
>>>>>>> 05d3630 (Optimized for reading only new data textfiles and attach them to the binary pickle.)
os.makedirs(dest_path + '30_plots/', exist_ok=True)

# plot SWE (Leica, Emlid, manual, laser, buoy, poles)
Expand All @@ -186,11 +144,7 @@

# plot SWE differences (Emlid, manual, laser, buoy, poles compared to Leica)
f.plot_all_diffSWE(dest_path, diffs_swe_daily, manual, laser_15min, buoy_daily, poles_daily,
<<<<<<< HEAD
save=save_plots, suffix='', leg=['Low-cost GNSS', 'Manual', 'Laser', 'Buoy', '_', '_', '_', '_', 'Poles'], y_lim=delta_swe_y_lim, x_lim=xlim_dates)
=======
save=save_plots, suffix='', leg=['Low-cost GNSS', 'Manual', 'Laser (SHM)'], y_lim=delta_swe_y_lim, x_lim=xlim_dates)
>>>>>>> 05d3630 (Optimized for reading only new data textfiles and attach them to the binary pickle.)

# plot boxplot of differences (Emlid, manual, laser compared to Leica)
f.plot_swediff_boxplot(dest_path, diffs_swe_daily, y_lim=delta_swe_y_lim, save=save_plots)
Expand All @@ -214,13 +168,8 @@
f.plot_SWE_density_acc(dest_path, swe_gnss_daily_leica.dropna(), swe_gnss_daily_emlid.dropna(), manual, laser_15min,
save=save_plots, std_leica=std_gnss_daily_leica.dropna(), std_emlid=std_gnss_daily_emlid.dropna(), suffix='', y_lim=acc_y_lim, x_lim=xlim_dates)

# Q: plot quality control
# plot number of satellites
<<<<<<< HEAD
f.plot_nrsat(dest_path, fil_df_leica.nr_sat, fil_df_emlid.nr_sat, save=save_plots, suffix='', y_lim=(0, 25), x_lim=xlim_dates)
=======
f.plot_nrsat(dest_path, fil_df_leica.nr_sat, fil_df_emlid.nr_sat, save=save_plots, suffix='', y_lim=(0, 35), x_lim=xlim_dates)
>>>>>>> 05d3630 (Optimized for reading only new data textfiles and attach them to the binary pickle.)

# plot ambiguity resolution state
f.plot_solquality(dest_path, df_enu_leica.amb_state, df_enu_emlid.amb_state, save=save_plots, suffix='', y_lim=(0, 100), x_lim=xlim_dates)
Expand All @@ -231,7 +180,6 @@


''' 8. Read GNSS-IR results '''
<<<<<<< HEAD
# read and filter gnss-ir snow accumulation results (processed using 'gnssrefl' on Linux)
df_rh = f.read_gnssir(dest_path, ubuntu_path, base_name, yy, copy=False, pickle='nmlb')
gnssir_acc, gnssir_acc_daily, gnssir_acc_daily_std, gnssir_rh_clean = f.filter_gnssir(df_rh, freq='2nd', threshold=2)
Expand Down Expand Up @@ -304,118 +252,9 @@
''' 13. Back up data '''
if solplot_backup is True:
# copy solution and plot directories back to server
=======
# best results: ele=5-30, f=2, azi=30-160 & 210-310

# read and filter gnss-ir snow accumulation results
df_rh, gnssir_acc, gnssir_acc_sel = f.read_gnssir(dest_path, ubuntu_path, base_name, yy, freq='2nd', excl_azi=True, copy=False, pickle='nmlb')

# plot gnss-ir snow accumulation results
gnssir_acc_median, gnssir_acc_std = f.plot_gnssir(dest_path, gnssir_acc_sel, laser_daily, leica_daily, emlid=None, manual=None, buoy=None, poles=None, leg=['GNSS-Reflectometry', '_', 'GNSS-Refractometry', 'Laser (SHM)'], save=save_plots, suffix='_leica', x_lim=xlim_dates)
f.plot_gnssir(dest_path, gnssir_acc_sel, laser_daily, leica_daily, emlid=emlid_daily, manual=None, buoy=None, poles=None, leg=['GNSS-Reflectometry', '_', 'High-end GNSS-Refractometry', 'Low-cost GNSS-Refractometry', 'Laser (SHM)'], save=save_plots, suffix='_emlid', x_lim=xlim_dates)


''' 9. Calculate snow density from GNSS-RR: Combine GNSS-IR & GNSS refractometry '''
density_leica, density_leica_cleaned = f.convert_swesh2density((leica_daily.dswe + 45).dropna(), (gnssir_acc_sel + 110).dropna())
density_emlid, density_emlid_cleaned = f.convert_swesh2density((emlid_daily.dswe + 45).dropna(), (gnssir_acc_sel + 110).dropna())

# calibrate with ref manual density above antenna where 1m is reached ['2022-09-22']
cal_leica_1m = manual.Density_aboveAnt['2022-09-22'] - density_leica_cleaned['2022-09-22']
cal_emlid_1m = manual.Density_aboveAnt['2022-09-22'] - density_emlid_cleaned['2022-09-22']

# plot density time series
f.plot_density(dest_path, density_leica_cleaned + cal_leica_1m, density_emlid=None, laser=None, manual=manual, leg=['High-end GNSS-RR', 'Manual'], save=save_plots, suffix='_leica', x_lim=xlim_dates)
f.plot_density(dest_path, density_leica_cleaned, density_emlid_cleaned, laser=None, manual=manual, leg=['High-end GNSS-RR', 'Low-cost GNSS-RR', 'Manual'], save=save_plots, suffix='_leica_emlid', x_lim=xlim_dates)
f.plot_density(dest_path, density_leica=None, density_emlid=density_emlid_cleaned + cal_emlid_1m, laser=None, manual=manual, leg=['_', 'Low-cost GNSS-RR', 'Manual'], save=save_plots, suffix='_emlid', x_lim=xlim_dates)

# plot difference to manual density observation
f.plot_diff_density(dest_path, manual.Density_aboveAnt - (density_leica_cleaned + cal_leica_1m), manual.Density_aboveAnt - (density_emlid_cleaned + cal_emlid_1m), laser=None, manual=manual, leg=['_', 'Low-cost GNSS-RR', 'Manual'], save=save_plots, suffix='_emlid', x_lim=xlim_dates)


''' 10. Back up data '''
if solplot_backup is True:
# copy solutions and plots directories back to server
>>>>>>> 05d3630 (Optimized for reading only new data textfiles and attach them to the binary pickle.)
f.copy_solplotsdirs(dest_path, scr_path + '../Processing/Run_RTKLib/data_neumayer/')

if total_backup is True:
# copy entire processing directory back to server
f.copy4backup(dest_path + '../', scr_path + '../Processing/Run_RTKLib/')
<<<<<<< HEAD


''' New stuff '''
# # STDs (AGM 2023)
# laser_15min.dsh_std.mean() # laser
# std_gnss_daily_leica.dropna().mean() # leica
# std_gnss_daily_emlid.dropna().mean() # emlid
# gnssir_acc.resample('D').std().median() # GNSS-IR
# rms
# np.sqrt(np.sum((f-g)**2)/n)
=======
>>>>>>> 05d3630 (Optimized for reading only new data textfiles and attach them to the binary pickle.)







''' New stuff '''
# TODO: test ppp diff - works!
# ppp_diff = ((df_ppp_rover.dh-df_ppp_ref.dh) * 1000)
# ppp_diff_sel = ppp_diff[(ppp_diff.index > '2021-12-01')]
# ppp_diff_sel = (ppp_diff_sel - ppp_diff_sel[0]).rolling('D').median()
#
# # Q: adjust for snow mast heightening (approx. 3m elevated several times a year)
# print('\ndata is corrected for snow mast heightening events (remove sudden jumps > 1m)')
# jump = ppp_diff_sel['2022-02-09'].median() # detect jumps (> 1000mm) in the dataset
#
# print('\njump of height %s is detected!' % jump)
# adj = ppp_diff_sel[(ppp_diff_sel.index >= '2022-02-09')] - jump # correct all observations after jump [0]
# ppp_diff_sel = pd.concat([ppp_diff_sel[~(ppp_diff_sel.index >= '2022-02-09')],
# adj]) # concatenate all original obs before jump with adjusted values after jump
#
# swe_gnss_ppp = ppp_diff_sel - ppp_diff_sel[0]
# swe_gnss_ppp = swe_gnss_ppp[(swe_gnss_ppp > -1000) & (swe_gnss_ppp < 1000)]
# swe_gnss_ppp_cleaned = swe_gnss_ppp[(swe_gnss_ppp.diff() > -50) & (swe_gnss_ppp.diff() < 50) ]
# swe_gnss_ppp_cleaned.plot();plt.grid(); plt.show()

# # STDs (AGM 2023)
# laser_15min.dsh_std.mean() # laser
# std_gnss_daily_leica.dropna().mean() # leica
# std_gnss_daily_emlid.dropna().mean() # emlid
# gnssir_acc.resample('D').std().median() # GNSS-IR
# rms
# np.sqrt(np.sum((f-g)**2)/n)

# TODO: add copy reflectometry solution files from ubuntu localhost
# Q: copy reflectometry solution files (*.txt) from the local Ubuntu server if not already existing
ubuntu_path = '//wsl.localhost/Ubuntu/home/sladina/test/gnssrefl/data/'
print(colored("\ncopy new reflectometry solution files", 'blue'))
print(ubuntu_path)
# get list of yearly directories newer than first year
for f in glob.glob(ubuntu_path + '2*'):
year = int(os.path.basename(f))
print(year)
# if year >= int('20' + yy):
# # copy missing laser observation files
# for f in glob.glob(laser_path + year + '/*.[ls]??'):
# file = os.path.basename(f)
# # skip files of 2021 before 26th nov (no gps data before installation)
# if int(file[2:8]) > 211125:
# if not os.path.exists(loc_laser_dir + file):
# shutil.copy2(f, loc_laser_dir)
# print("file copied from %s to %s" % (f, loc_laser_dir))
# else:
# # print(colored("\nfile in destination already exists: %s, \ncopy aborted!!!" % dest_path, 'yellow'))
# pass
# else:
# pass
# else:
# pass
print(colored("\nnew reflectometry solution files copied", 'blue'))

# test
if jump2.iloc[0] is not None:
print(jump2.iloc[0], jump2.index.format()[0])

0 comments on commit d4f67ac

Please sign in to comment.