Skip to content

Commit

Permalink
Merge pull request #36 from paarnes/dev
Browse files Browse the repository at this point in the history
Dev
  • Loading branch information
paarnes authored Mar 17, 2024
2 parents 8bb14a9 + 1ca51b4 commit 53b58e2
Show file tree
Hide file tree
Showing 6 changed files with 232 additions and 38 deletions.
Binary file added Results_example/GPS_results.xlsx
Binary file not shown.
Binary file added Results_example/result_table.PNG
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
98 changes: 60 additions & 38 deletions src/gnssmultipath/GNSS_MultipathAnalysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,11 +13,13 @@
import numpy as np
from tqdm import tqdm
from gnssmultipath.readRinexObs import readRinexObs
from gnssmultipath.Geodetic_functions import gpstime_to_utc_datefmt
from gnssmultipath.computeSatElevations import computeSatElevations
from gnssmultipath.computeSatElevAzimuth_fromNav import computeSatElevAzimuth_fromNav
from gnssmultipath.signalAnalysis import signalAnalysis
from gnssmultipath.detectClockJumps import detectClockJumps
from gnssmultipath.writeOutputFile import writeOutputFile
from gnssmultipath.createCSVfile import createCSVfile
from gnssmultipath.make_polarplot import make_polarplot,make_skyplot, make_polarplot_SNR, plot_SNR_wrt_elev
from gnssmultipath.make_polarplot_dont_use_TEX import make_polarplot_dont_use_TEX, make_skyplot_dont_use_TEX, make_polarplot_SNR_dont_use_TEX, plot_SNR_wrt_elev_dont_use_TEX
from gnssmultipath.plotResults import plotResults, plotResults_dont_use_TEX, make_barplot, make_barplot_dont_use_TEX
Expand Down Expand Up @@ -45,6 +47,8 @@ def GNSS_MultipathAnalysis(rinObsFilename: str,
include_SNR: bool = True,
save_results_as_pickle: bool = True,
save_results_as_compressed_pickle: bool = False,
write_results_to_csv: bool = True,
output_csv_delimiter: str = ';',
nav_data_rate: int = 60,
includeResultSummary: Union[bool, None] = None,
includeCompactSummary: Union[bool, None] = None,
Expand Down Expand Up @@ -134,6 +138,10 @@ def GNSS_MultipathAnalysis(rinObsFilename: str,
save_results_as_compressed_pickle : boolean. If True, the results will be stored as dictionary in form of a binary compressed pickle file (zstd compression). Default set to False.
write_results_to_csv: boolean. If True, a subset of the results will be exported as a CSV file. Default is True.
output_csv_delimiter: str. Set the delimiter of the CSV file. Default is semi colon (;).
nav_data_rate: integer. The desired data rate of ephemerides given in minutes. Default is 60 min. The purpose with this
Expand Down Expand Up @@ -630,17 +638,19 @@ def GNSS_MultipathAnalysis(rinObsFilename: str,
rinex_obs_filename = rinObsFilename.split('/')
rinex_obs_filename = rinex_obs_filename[-1]
analysisResults['ExtraOutputInfo'] = {}
analysisResults['ExtraOutputInfo']['rinex_obs_filename'] = rinex_obs_filename
analysisResults['ExtraOutputInfo']['markerName'] = markerName
analysisResults['ExtraOutputInfo']['rinexVersion'] = rinexVersion
analysisResults['ExtraOutputInfo']['rinexProgr'] = rinexProgr
analysisResults['ExtraOutputInfo']['recType'] = recType
analysisResults['ExtraOutputInfo']['tFirstObs'] = tFirstObs
analysisResults['ExtraOutputInfo']['tLastObs'] = tLastObs
analysisResults['ExtraOutputInfo']['tInterval'] = tInterval
analysisResults['ExtraOutputInfo']['GLO_Slot2ChannelMap'] = GLO_Slot2ChannelMap
analysisResults['ExtraOutputInfo']['nEpochs'] = nepochs
analysisResults['ExtraOutputInfo']['elevation_cutoff'] = cutoff_elevation_angle
analysisResults['ExtraOutputInfo']['rinex_obs_filename'] = rinex_obs_filename
analysisResults['ExtraOutputInfo']['markerName'] = markerName
analysisResults['ExtraOutputInfo']['rinexVersion'] = rinexVersion
analysisResults['ExtraOutputInfo']['rinexProgr'] = rinexProgr
analysisResults['ExtraOutputInfo']['recType'] = recType
analysisResults['ExtraOutputInfo']['tFirstObs'] = tFirstObs
analysisResults['ExtraOutputInfo']['tLastObs'] = tLastObs
analysisResults['ExtraOutputInfo']['tInterval'] = tInterval
analysisResults['ExtraOutputInfo']['time_epochs_gps_time'] = time_epochs
analysisResults['ExtraOutputInfo']['time_epochs_utc_time'] = gpstime_to_utc_datefmt(time_epochs)
analysisResults['ExtraOutputInfo']['GLO_Slot2ChannelMap'] = GLO_Slot2ChannelMap
analysisResults['ExtraOutputInfo']['nEpochs'] = nepochs
analysisResults['ExtraOutputInfo']['elevation_cutoff'] = cutoff_elevation_angle

## -- Store default limits or user set limits in dict
if phaseCodeLimit == 0:
Expand Down Expand Up @@ -679,6 +689,8 @@ def GNSS_MultipathAnalysis(rinObsFilename: str,
outputFilename = baseFileName.split('.')[0] + '_Report.txt'
writeOutputFile(outputFilename, outputDir, analysisResults, includeResultSummary, includeCompactSummary, includeObservationOverview, includeLLIOverview)
print('INFO: The output file %s has been written.\n' % (outputFilename))


## -- Make barplot if plotEstimates is True
if plotEstimates:
print('INFO: Making bar plot. Please wait...\n')
Expand Down Expand Up @@ -721,48 +733,58 @@ def GNSS_MultipathAnalysis(rinObsFilename: str,
else:
make_polarplot_dont_use_TEX(analysisResults, graphDir)

if include_SNR:
# Seaching for SNR codes
for sys in GNSS_obs.keys():
GNSSsystemIndex = [k for k in GNSSsystems if GNSSsystems[k] == sys][0]
SNR_codes = [SNR_code for SNR_code in obsCodes[GNSSsystemIndex][sys] if 'S' in SNR_code[0]]
if len(SNR_codes) != 0:
print('INFO: Making a plot of the Signal To Noise Ration (SNR). Please wait ...')
if use_LaTex:
try:
make_polarplot_SNR(analysisResults,GNSS_obs,GNSSsystems,obsCodes,graphDir)
plot_SNR_wrt_elev(analysisResults,GNSS_obs,GNSSsystems,obsCodes,graphDir,tInterval)
except:
make_polarplot_SNR_dont_use_TEX(analysisResults,GNSS_obs,GNSSsystems,obsCodes,graphDir)
plot_SNR_wrt_elev_dont_use_TEX(analysisResults,GNSS_obs,GNSSsystems,obsCodes,graphDir,tInterval)
else:

if include_SNR:
# Seaching for SNR codes
for sys in GNSS_obs.keys():
GNSSsystemIndex = [k for k in GNSSsystems if GNSSsystems[k] == sys][0]
SNR_codes = [SNR_code for SNR_code in obsCodes[GNSSsystemIndex][sys] if 'S' in SNR_code[0]]
if plotEstimates and len(SNR_codes) != 0:
print('INFO: Making a plot of the Signal To Noise Ration (SNR). Please wait ...')
if use_LaTex:
try:
make_polarplot_SNR(analysisResults,GNSS_obs,GNSSsystems,obsCodes,graphDir)
plot_SNR_wrt_elev(analysisResults,GNSS_obs,GNSSsystems,obsCodes,graphDir,tInterval)
except:
make_polarplot_SNR_dont_use_TEX(analysisResults,GNSS_obs,GNSSsystems,obsCodes,graphDir)
plot_SNR_wrt_elev_dont_use_TEX(analysisResults,GNSS_obs,GNSSsystems,obsCodes,graphDir,tInterval)
else:
logger.warning("INFO(GNSS_MultipathAnalysis): There is no SNR codes available in the RINEX files. Plot of the Signal To Noise Ration is not possible.")
make_polarplot_SNR_dont_use_TEX(analysisResults,GNSS_obs,GNSSsystems,obsCodes,graphDir)
plot_SNR_wrt_elev_dont_use_TEX(analysisResults,GNSS_obs,GNSSsystems,obsCodes,graphDir,tInterval)
else:
logger.warning("INFO(GNSS_MultipathAnalysis): There is no SNR codes available in the RINEX files. Plot of the Signal To Noise Ration is not possible.")

for syst in GNSSsystems.keys():
curr_syst =GNSSsystemCode2Fullname[GNSSsystems[syst]]
analysisResults[curr_syst]['SNR'] = {}
curr_obscodes = obsCodes[syst][GNSSsystems[syst]]
snr_codes_idx = [n for n, l in enumerate(curr_obscodes) if l.startswith('S')]
for code_idx in snr_codes_idx:
signal = curr_obscodes[code_idx]
curr_ban = [element for element in analysisResults[curr_syst].keys() if element.endswith(signal[1])][0]
curr_signal = np.stack(list(GNSS_obs[GNSSsystems[syst]].values()))[:, :, code_idx]
curr_signal = np.squeeze(curr_signal)
analysisResults[curr_syst]['SNR'][signal] = curr_signal

for syst in GNSSsystems.keys():
curr_syst =GNSSsystemCode2Fullname[GNSSsystems[syst]]
analysisResults[curr_syst]['SNR'] = {}
curr_obscodes = obsCodes[syst][GNSSsystems[syst]]
snr_codes_idx = [n for n, l in enumerate(curr_obscodes) if l.startswith('S')]
for code_idx in snr_codes_idx:
signal = curr_obscodes[code_idx]
curr_ban = [element for element in analysisResults[curr_syst].keys() if element.endswith(signal[1])][0]
curr_signal = np.stack(list(GNSS_obs[GNSSsystems[syst]].values()))[:, :, code_idx]
curr_signal = np.squeeze(curr_signal)
analysisResults[curr_syst]['SNR'][signal] = curr_signal

if write_results_to_csv:
result_dir = os.path.join(outputDir, "Result_files_CSV")
if not os.path.exists(result_dir):
os.makedirs(result_dir)
createCSV = createCSVfile(analysisResults, result_dir, output_csv_delimiter)
createCSV.write_results_to_csv()


## -- Saving the workspace as a binary pickle file ---
if save_results_as_pickle:
pickle_filename = 'analysisResults.pkl'
print(f'\nINFO: The analysis results are being written to the file {pickle_filename}. Please wait..')
results_name = os.path.join(outputDir, pickle_filename)
PickleHandler.write_zstd_pickle(analysisResults, results_name)
print(f'INFO: The analysis results has been written to the file {pickle_filename}.\n')
elif save_results_as_compressed_pickle:
pickle_filename = 'analysisResults.pkl.zst'
print(f'\nINFO: The analysis results are being written to the file {pickle_filename}. Please wait..')
results_name = os.path.join(outputDir, pickle_filename)
PickleHandler.write_zstd_pickle(analysisResults, results_name)
print(f'INFO: The analysis results has been written to the file {pickle_filename}.\n')
Expand Down
21 changes: 21 additions & 0 deletions src/gnssmultipath/Geodetic_functions.py
Original file line number Diff line number Diff line change
Expand Up @@ -632,6 +632,27 @@ def gpstime2date_arrays(week: Union[List[int], np.ndarray], tow: Union[List[floa
return date_array


def convert_to_datetime_vectorized(time_datetime: np.ndarray) -> list:
"""
Convert numpy array of datetime components to datetime strings with specified format.
"""
# Convert datetime components to datetime objects
datetimes = [datetime(*components) for components in time_datetime]
# Format datetime objects to strings
datetime_strings = np.array([dt.strftime("%Y-%m-%d %H:%M:%S") for dt in datetimes]).tolist()
return datetime_strings


def gpstime_to_utc_datefmt(time_epochs_gpstime: np.ndarray) -> list:
"""
Coverters form GPS time to UTC with formatting.
Ex output: "2022-01-01 02:23:30".
"""
time_datetime = gpstime2date_arrays(*time_epochs_gpstime.T)
return convert_to_datetime_vectorized(time_datetime)



def date2gpstime_vectorized(gregorian_date_array):
"""
Expand Down
150 changes: 150 additions & 0 deletions src/gnssmultipath/createCSVfile.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,150 @@
import pandas as pd
import numpy as np
import os
import sys



class createCSVfile:
def __init__(self, analysisResults, output_dir, column_delimiter=';'):
self.analysisResults = analysisResults
self.output_dir = output_dir
self.column_delimiter = column_delimiter
self.results_dict = {}
self.format_rules = {
"Azimuth": "{:.2f}",
"Elevation": "{:.2f}",
"MP_": "{:.4f}",
"SNR_": "{:.1f}"
}
self.time_stamps = self.analysisResults["ExtraOutputInfo"]["time_epochs_utc_time"]
self.mp_data_lst = ["PRN; Time_UTC; Elevation; Azimuth"]
self.GNSSsystemCode2Fullname = {'G': 'GPS', 'R': 'GLONASS', 'E': 'Galileo', 'C': 'BeiDou'}
self.GNSS_Name2Code = {v: k for k, v in self.GNSSsystemCode2Fullname.items()}
self.results_dict = self.build_results_dict()


def flatten_result_array(self, arr):
"""
Flatten a numpy array to 1D
"""
flatten_array = arr[:, 1:].T.ravel().tolist()
return flatten_array

def set_float_fmt_dataframe(self, df):
for column in df.columns:
for prefix, fmt in self.format_rules.items():
if column.startswith(prefix):
df[column] = df[column].map(lambda x: fmt.format(x))
break

def extract_multipath_and_put_in_result_dict(self, sys_name):
"""
Extract the multipath values from "analysisResults" dictionary
and gather them in the results_dict.
"""

for band in self.analysisResults[sys_name]["Bands"]:
for code in self.analysisResults[sys_name][band]["Codes"]:
curr_code_data = self.analysisResults[sys_name][band].get(code, None) if not isinstance(code, list) else None
if curr_code_data is not None:
rms_multipath_avg = np.round(curr_code_data["multipath_range1"], 4)
else:
continue
mp_header = f"MP_{code}"
self.mp_data_lst.append(mp_header)
self.results_dict[sys_name][mp_header] = rms_multipath_avg


def extract_SNR_and_put_in_result_dict(self, sys_name):
"""
Extract the SNR values from "analysisResults" dictionary
and gather them in the results_dict.
"""
SNR_dict = self.analysisResults[sys_name].get("SNR", None)
if SNR_dict is not None:
for signal_code, snr_array in SNR_dict.items():
snr_array[snr_array == 0] = np.nan # convert null to np.nan
if not np.all(np.isnan(snr_array)):
signal_header = f"SNR_{signal_code}"
self.mp_data_lst.append(signal_header)
# Append SNR data to results_dict if not all elements are np.nan
self.results_dict[sys_name][signal_header] = snr_array



def build_results_dict(self):
# results_dict = {}
GNSS_systems = list(self.analysisResults["Sat_position"].keys())

for gnss_sys in GNSS_systems:
sys_name = self.GNSSsystemCode2Fullname[gnss_sys]
self.results_dict[sys_name] = {"Elevation": [], "Azimuth": []}

# Build up result dict
for gnss_sys in GNSS_systems:
curr_sys = self.analysisResults["Sat_position"][gnss_sys]
sys_name = self.GNSSsystemCode2Fullname[gnss_sys]
self.results_dict[sys_name]["Azimuth"] = curr_sys["azimuth"]
self.results_dict[sys_name]["Elevation"] = curr_sys["elevation"]

# Extract multipath and SNR values and store in results_dict
self.extract_multipath_and_put_in_result_dict(sys_name)
self.extract_SNR_and_put_in_result_dict(sys_name)


return self.results_dict

def write_results_to_csv(self):
for sys_name, sys_data in self.results_dict.items():
# Extract data into arrays
sys_code = self.GNSS_Name2Code[sys_name]
timestamps = self.time_stamps * sys_data["Elevation"][:, 1:].shape[1]

prns = list(range(1, sys_data["Azimuth"].shape[1]))
prns = [f"{sys_code}{prn:02d}" for prn in prns]
prn_repeated = list(np.repeat(prns, sys_data["Azimuth"].shape[0]))

# Flatten numpy array to 1D
az = self.flatten_result_array(np.round(sys_data["Azimuth"], 2))
el = self.flatten_result_array(np.round(sys_data["Elevation"], 2))

# Create a DataFrame for the current system
df = pd.DataFrame({
"PRN": prn_repeated,
"Time_UTC": timestamps,
"Azimuth": az,
"Elevation": el
}, dtype=object)


# Add SNR columns to the DataFrame
if any(key.startswith("MP") for key in sys_data.keys()):
mp_headers = [header for header in sys_data.keys() if header.startswith("MP_")]
for i, header in enumerate(mp_headers):
df[header] = self.flatten_result_array(sys_data[header])

# Add SNR columns to the DataFrame
if any(key.startswith("SNR") for key in sys_data.keys()):
snr_headers = [header for header in sys_data.keys() if header.startswith("SNR_")]
for i, header in enumerate(snr_headers):
df[header] = self.flatten_result_array(sys_data[header])

# Set specified float formatting on the columns
self.set_float_fmt_dataframe(df)

# Remove rows where all values except PRN and time are np.nan
df[df.columns[2:]] = df[df.columns[2:]].apply(pd.to_numeric, errors='coerce') # Convert selected columns to numeric
df.dropna(subset=['Elevation'], inplace=True) # drop rows where the satellite is below the horizon
output_file = os.path.join(self.output_dir, f"{sys_name}_results.csv")
print(f'INFO: The result CSV file {output_file} has been written')
df.to_csv(output_file, index=False, sep=self.column_delimiter)


if __name__ =="__main__":
sys.path.append(r"C:\Users\perhe\OneDrive\Documents\Python_skript\GNSS_repo\src")
from gnssmultipath import PickleHandler
analysisResults = PickleHandler.read_zstd_pickle(r"C:\Users\perhe\Desktop\CSV_export\analysisResults.pkl")
outputDir = r"C:\Users\perhe\Desktop\CSV_export\TEST"
createCSV = createCSVfile(analysisResults, outputDir)
createCSV.write_results_to_csv()
1 change: 1 addition & 0 deletions tests/test_GNSS_MultipathAnalysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,7 @@ def test_GNSS_MultipathAnalysis_broadcast_navfile():
result = GNSS_MultipathAnalysis(rinObsFilename=rinObs_file, broadcastNav1=broadNav_file,
plotEstimates=False,
plot_polarplot=False,
write_results_to_csv=False,
nav_data_rate=120)

expected_result = PickleHandler.read_zstd_pickle(expected_res)
Expand Down

0 comments on commit 53b58e2

Please sign in to comment.