From 47878166a79a0b5bbf61c35edfd7993f3c9a3235 Mon Sep 17 00:00:00 2001 From: adeane-ga Date: Thu, 5 Jan 2023 12:19:14 +1100 Subject: [PATCH 1/6] Add function to remove velocity from sinex --- geodepy/gnss.py | 160 +++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 159 insertions(+), 1 deletion(-) diff --git a/geodepy/gnss.py b/geodepy/gnss.py index e5b510e..a8aa694 100644 --- a/geodepy/gnss.py +++ b/geodepy/gnss.py @@ -8,8 +8,9 @@ """ from datetime import datetime -from numpy import zeros +from numpy import zeros, delete from geodepy.angles import DMSAngle +import re def list_sinex_blocks(file): @@ -699,3 +700,160 @@ def remove_stns_sinex(sinex, sites): out.write('%ENDSNX\n') return + +def remove_velocity_sinex(sinex): + """This function reads in a SINEX file and removes the + velocity parameters, including the zeros of the VCV matrix, + + :param str sinex: input SINEX file + return: SINEX file output.snx + """ + + # From header, confirm that the SINEX has velocity parameters + header = read_sinex_header_line(sinex) + if header[70:71] == 'V': + pass + else: + print("Not removing velocities because SINEX file does not have velocity parameters.") + exit() + + # Open the output file + with open('output.snx', 'w') as out: + + # With header line: + # - update the creation time + # - update number of parameter estimates + # - then write to file + old_creation_time = header[15:27] + creation_time = set_creation_time() + header = header.replace(old_creation_time, creation_time) + old_num_params = header[60:65] + num_params = int(old_num_params) / 2 + header = header.replace(old_num_params, str(num_params)) + out.write(header) + del header + + # Read in the +SITE/ID block and write to file + site_id = read_sinex_site_id_block(sinex) + for line in site_id: + out.write(f"{line}") + del site_id + + # Read in the +SOLUTION/EPOCHS block and write to file + # - also collect count on number of sites for use later + numSites = 0 + solution_epochs = read_sinex_solution_epochs_block(sinex) + for line in solution_epochs: + out.write(f"{line}") + if line[0]!="+" and line[0]!="*" and line[0]!="-": + numSites+=1 + del solution_epochs + + # Read in the +SOLUTION/ESTIMATE block: + # - gather velocity indices + # - remove velocity rows + # - change indices to account for removed velocities + # - write to file + vel_indices = [] + estimate_number = 0 + solution_estimate = read_sinex_solution_estimate_block(sinex) + for line in solution_estimate: + if line[7:10]=="VEL": + vel_indices.append(int(line[0:6])) + elif line[0]=="+" or line[0]=="*" or line[0]=="-": + out.write(f"{line}") + else: + estimate_number+=1 + number = '{:5d}'.format(estimate_number) + line = ' ' + number + line[6:] + out.write(f"{line}") + del solution_estimate + + # Read in the +SOLUTION/MATRIX_ESTIMATE block + # - identify if matrix is upper or lower triangle form + # - form full VCV matrix + # - remove all rows/columns associated with velocity + # - write to file with new matrix indices (para1, para2) + solution_matrix_estimate = read_sinex_solution_matrix_estimate_block(sinex) + if solution_matrix_estimate[0][26:27] == 'L': + matrix = 'lower' + elif solution_matrix_estimate[0][26:27] == 'U': + matrix = 'upper' + # Size of original VCV matrix + matrix_dim = numSites*6 + # Setup matrix of zeros + Q = zeros((matrix_dim, matrix_dim)) + # Write initial comment line(s), save last comment line, and form VCV matrix + for line in solution_matrix_estimate: + if line[0]=="+": + out.write(f"{line}") + continue + if line[0]=="*": + out.write(f"{line}") + continue + if line[0]=="-": + block_end_title = line + else: + lineVCV = re.split('\s+', line) + lineVCV = list(filter(None, lineVCV)) + row = int(lineVCV[0]) + col = int(lineVCV[1]) + if len(lineVCV) >= 3: + q1 = float(lineVCV[2]) + Q[row-1, col-1] = q1 + Q[col-1, row-1] = q1 + if len(lineVCV) >= 4: + q2 = float(lineVCV[3]) + Q[row-1, col] = q2 + Q[col, row-1] = q2 + if len(lineVCV) >= 5: + q3 = float(lineVCV[4]) + Q[row-1, col+1] = q3 + Q[col+1, row-1] = q3 + # Remove velocity rows/columns from VCV matrix + num_removed = 0 + for i in vel_indices: + Q = delete(Q, i-1-num_removed, 0) + Q = delete(Q, i-1-num_removed, 1) + num_removed += 1 + # Write VCV matrix to SINEX (lower triangle) + if matrix == 'lower': + i = 0 + j = 0 + while i < len(Q): + while j <= i: + out.write(f" {i+1:5d} {j+1:5d} {Q[i,j]:21.14E} ") + j += 1 + if j <= i: + out.write(f"{Q[i,j]:21.14E} ") + j += 1 + if j <= i: + out.write(f"{Q[i,j]:21.14E}") + j += 1 + out.write(" \n") + j = 0 + i += 1 + # Write VCV matrix to SINEX (upper triangle) + if matrix == 'upper': + j = 0 + for i in range(len(Q)): + j = i + while j < len(Q): + out.write(f" {i+1:5d} {j+1:5d} {Q[i,j]:21.14E} ") + j += 1 + if j < len(Q): + out.write(f"{Q[i,j]:21.14E} ") + j += 1 + if j < len(Q): + out.write(f"{Q[i,j]:21.14E}") + j += 1 + out.write(" \n") + # Write out end of block line, and delete large variables + out.write(block_end_title) + del solution_matrix_estimate + del Q + + # Write out the trailer line + out.write('%ENDSNX\n') + + return From 403e008806189d0cce0cac5c6995056851bcee21 Mon Sep 17 00:00:00 2001 From: adeane-ga Date: Thu, 5 Jan 2023 15:16:59 +1100 Subject: [PATCH 2/6] Minor refinement of remove_velocity_sinex() --- geodepy/gnss.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/geodepy/gnss.py b/geodepy/gnss.py index a8aa694..cdff527 100644 --- a/geodepy/gnss.py +++ b/geodepy/gnss.py @@ -722,7 +722,8 @@ def remove_velocity_sinex(sinex): # With header line: # - update the creation time - # - update number of parameter estimates + # - update number of parameter estimates + # - remove 'V' from parameter list # - then write to file old_creation_time = header[15:27] creation_time = set_creation_time() @@ -730,6 +731,7 @@ def remove_velocity_sinex(sinex): old_num_params = header[60:65] num_params = int(old_num_params) / 2 header = header.replace(old_num_params, str(num_params)) + header = header.replace('V', '') out.write(header) del header @@ -769,7 +771,7 @@ def remove_velocity_sinex(sinex): out.write(f"{line}") del solution_estimate - # Read in the +SOLUTION/MATRIX_ESTIMATE block + # Read in the +SOLUTION/MATRIX_ESTIMATE block: # - identify if matrix is upper or lower triangle form # - form full VCV matrix # - remove all rows/columns associated with velocity @@ -792,7 +794,7 @@ def remove_velocity_sinex(sinex): out.write(f"{line}") continue if line[0]=="-": - block_end_title = line + block_end = line else: lineVCV = re.split('\s+', line) lineVCV = list(filter(None, lineVCV)) @@ -849,7 +851,7 @@ def remove_velocity_sinex(sinex): j += 1 out.write(" \n") # Write out end of block line, and delete large variables - out.write(block_end_title) + out.write(block_end) del solution_matrix_estimate del Q From ac91b9b2c55ab23877e2f317a8c0bb4d324a8a4b Mon Sep 17 00:00:00 2001 From: adeane-ga Date: Thu, 5 Jan 2023 15:36:05 +1100 Subject: [PATCH 3/6] Add function to only remove zeros from matrix --- geodepy/gnss.py | 99 ++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 82 insertions(+), 17 deletions(-) diff --git a/geodepy/gnss.py b/geodepy/gnss.py index cdff527..527c43c 100644 --- a/geodepy/gnss.py +++ b/geodepy/gnss.py @@ -703,7 +703,7 @@ def remove_stns_sinex(sinex, sites): def remove_velocity_sinex(sinex): """This function reads in a SINEX file and removes the - velocity parameters, including the zeros of the VCV matrix, + velocity parameters, including the zeros of the SOLUTION/MATRIX_ESTIMATE block, :param str sinex: input SINEX file return: SINEX file output.snx @@ -773,7 +773,7 @@ def remove_velocity_sinex(sinex): # Read in the +SOLUTION/MATRIX_ESTIMATE block: # - identify if matrix is upper or lower triangle form - # - form full VCV matrix + # - form full matrix # - remove all rows/columns associated with velocity # - write to file with new matrix indices (para1, para2) solution_matrix_estimate = read_sinex_solution_matrix_estimate_block(sinex) @@ -781,11 +781,11 @@ def remove_velocity_sinex(sinex): matrix = 'lower' elif solution_matrix_estimate[0][26:27] == 'U': matrix = 'upper' - # Size of original VCV matrix + # Size of original matrix matrix_dim = numSites*6 # Setup matrix of zeros Q = zeros((matrix_dim, matrix_dim)) - # Write initial comment line(s), save last comment line, and form VCV matrix + # Write initial comment line(s), save last comment line, and form matrix for line in solution_matrix_estimate: if line[0]=="+": out.write(f"{line}") @@ -796,29 +796,29 @@ def remove_velocity_sinex(sinex): if line[0]=="-": block_end = line else: - lineVCV = re.split('\s+', line) - lineVCV = list(filter(None, lineVCV)) - row = int(lineVCV[0]) - col = int(lineVCV[1]) - if len(lineVCV) >= 3: - q1 = float(lineVCV[2]) + lineMAT = re.split('\s+', line) + lineMAT = list(filter(None, lineMAT)) + row = int(lineMAT[0]) + col = int(lineMAT[1]) + if len(lineMAT) >= 3: + q1 = float(lineMAT[2]) Q[row-1, col-1] = q1 Q[col-1, row-1] = q1 - if len(lineVCV) >= 4: - q2 = float(lineVCV[3]) + if len(lineMAT) >= 4: + q2 = float(lineMAT[3]) Q[row-1, col] = q2 Q[col, row-1] = q2 - if len(lineVCV) >= 5: - q3 = float(lineVCV[4]) + if len(lineMAT) >= 5: + q3 = float(lineMAT[4]) Q[row-1, col+1] = q3 Q[col+1, row-1] = q3 - # Remove velocity rows/columns from VCV matrix + # Remove velocity rows/columns from matrix num_removed = 0 for i in vel_indices: Q = delete(Q, i-1-num_removed, 0) Q = delete(Q, i-1-num_removed, 1) num_removed += 1 - # Write VCV matrix to SINEX (lower triangle) + # Write matrix to SINEX (lower triangle) if matrix == 'lower': i = 0 j = 0 @@ -835,7 +835,7 @@ def remove_velocity_sinex(sinex): out.write(" \n") j = 0 i += 1 - # Write VCV matrix to SINEX (upper triangle) + # Write matrix to SINEX (upper triangle) if matrix == 'upper': j = 0 for i in range(len(Q)): @@ -859,3 +859,68 @@ def remove_velocity_sinex(sinex): out.write('%ENDSNX\n') return + +def remove_matrixzeros_sinex(sinex): + """This function reads in a SINEX file and removes the + zeros from the SOLUTION/MATRIX_ESTIMATE block only. + + :param str sinex: input SINEX file + return: SINEX file output.snx + """ + + sinex = "apref20220423_NT.snx" + #sinex = "apref20210410_LowerTriangle.snx" + + # Open the output file + with open('output_NT.snx', 'w') as out: + + # With header line: + # - update the creation time + # - then write to file + header = read_sinex_header_line(sinex) + old_creation_time = header[15:27] + creation_time = set_creation_time() + header = header.replace(old_creation_time, creation_time) + out.write(header) + del header + + # Read in the +SITE/ID block and write to file + site_id = read_sinex_site_id_block(sinex) + for line in site_id: + out.write(f"{line}") + del site_id + + # Read in the +SOLUTION/EPOCHS block and write to file + solution_epochs = read_sinex_solution_epochs_block(sinex) + for line in solution_epochs: + out.write(f"{line}") + del solution_epochs + + # Read in the +SOLUTION/ESTIMATE block + solution_estimate = read_sinex_solution_estimate_block(sinex) + for line in solution_estimate: + out.write(f"{line}") + del solution_estimate + + # Read in the +SOLUTION/MATRIX_ESTIMATE block: + # - Remove lines that contain only zeros + solution_matrix_estimate = read_sinex_solution_matrix_estimate_block(sinex) + for line in solution_matrix_estimate: + col = line.split() + numCol = len(col) + if numCol==3: + if col[2]=="0.00000000000000e+00": + continue + if numCol==4: + if col[2]=="0.00000000000000e+00" and col[3]=="0.00000000000000e+00": + continue + if numCol==5: + if col[2]=="0.00000000000000e+00" and col[3]=="0.00000000000000e+00" and col[4]=="0.00000000000000e+00": + continue + out.write(line) + del solution_matrix_estimate + + # Write out the trailer line + out.write('%ENDSNX\n') + + return From 495629a8f2f01714a60e3b771a0b0a533576a5f4 Mon Sep 17 00:00:00 2001 From: adeane-ga Date: Thu, 5 Jan 2023 15:50:05 +1100 Subject: [PATCH 4/6] Minor edit to remove_matrixzeros_sinex() --- geodepy/gnss.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/geodepy/gnss.py b/geodepy/gnss.py index 527c43c..88a22b6 100644 --- a/geodepy/gnss.py +++ b/geodepy/gnss.py @@ -867,12 +867,9 @@ def remove_matrixzeros_sinex(sinex): :param str sinex: input SINEX file return: SINEX file output.snx """ - - sinex = "apref20220423_NT.snx" - #sinex = "apref20210410_LowerTriangle.snx" - + # Open the output file - with open('output_NT.snx', 'w') as out: + with open('output.snx', 'w') as out: # With header line: # - update the creation time From ee5c05725f4451182fdbf901010ab90115b76c94 Mon Sep 17 00:00:00 2001 From: adeane-ga Date: Fri, 6 Jan 2023 11:34:55 +1100 Subject: [PATCH 5/6] Refine SINEX formatting, and save SINEX comments --- geodepy/gnss.py | 67 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) diff --git a/geodepy/gnss.py b/geodepy/gnss.py index 88a22b6..9a94883 100644 --- a/geodepy/gnss.py +++ b/geodepy/gnss.py @@ -44,6 +44,28 @@ def print_sinex_comments(file): if line.startswith('-FILE/COMMENT'): go = False +def read_sinex_comments(file): + """This function reads comments in a SINEX file. + + :param str file: the input SINEX file + :return: comments block + :rtype: list of strings + """ + comments = [] + go = False + with open(file) as f: + for line in f.readlines(): + if line.startswith('+FILE/COMMENT'): + go = True + if go: + comments.append(line.strip()) + if line.startswith('-FILE/COMMENT'): + go = False + + comments.insert(-1,f"* File created by Geodepy.gnss.py at {datetime.now().strftime('%d-%m-%Y, %H:%M')}") + + return comments + def set_creation_time(): """This function sets the creation time, in format YY:DDD:SSSSS, for use @@ -588,6 +610,15 @@ def remove_stns_sinex(sinex, sites): header = header.replace(str(old_num_params), str(num_params)) out.write(header) + out.write("*-------------------------------------------------------------------------------\n") + + # Read in the +FILE/COMMENT block and write to output file + comments = read_sinex_comments(sinex) + for i in comments: + out.write(f"{i}\n") + + out.write("*-------------------------------------------------------------------------------\n") + # Read in the site ID block and write out the sites not being removed site_id = read_sinex_site_id_block(sinex) for line in site_id: @@ -600,6 +631,8 @@ def remove_stns_sinex(sinex, sites): out.write(line) del site_id + out.write("*-------------------------------------------------------------------------------\n") + # Read in the solution epochs block and write out the epochs of the # sites not being removed solution_epochs = read_sinex_solution_epochs_block(sinex) @@ -613,6 +646,8 @@ def remove_stns_sinex(sinex, sites): out.write(line) del solution_epochs + out.write("*-------------------------------------------------------------------------------\n") + # Read in the solution estimate block and write out the estimates of # the sites not being removed skip = [] @@ -634,6 +669,8 @@ def remove_stns_sinex(sinex, sites): out.write(line) del solution_estimate + out.write("*-------------------------------------------------------------------------------\n") + # Read in the matrix estimate block and write out minus the sites # being removed vcv = {} @@ -734,12 +771,23 @@ def remove_velocity_sinex(sinex): header = header.replace('V', '') out.write(header) del header + + out.write("*-------------------------------------------------------------------------------\n") + + # Read in the +FILE/COMMENT block and write to output file + comments = read_sinex_comments(sinex) + for i in comments: + out.write(f"{i}\n") + + out.write("*-------------------------------------------------------------------------------\n") # Read in the +SITE/ID block and write to file site_id = read_sinex_site_id_block(sinex) for line in site_id: out.write(f"{line}") del site_id + + out.write("*-------------------------------------------------------------------------------\n") # Read in the +SOLUTION/EPOCHS block and write to file # - also collect count on number of sites for use later @@ -751,6 +799,8 @@ def remove_velocity_sinex(sinex): numSites+=1 del solution_epochs + out.write("*-------------------------------------------------------------------------------\n") + # Read in the +SOLUTION/ESTIMATE block: # - gather velocity indices # - remove velocity rows @@ -771,6 +821,8 @@ def remove_velocity_sinex(sinex): out.write(f"{line}") del solution_estimate + out.write("*-------------------------------------------------------------------------------\n") + # Read in the +SOLUTION/MATRIX_ESTIMATE block: # - identify if matrix is upper or lower triangle form # - form full matrix @@ -880,12 +932,23 @@ def remove_matrixzeros_sinex(sinex): header = header.replace(old_creation_time, creation_time) out.write(header) del header + + out.write("*-------------------------------------------------------------------------------\n") + + # Read in the +FILE/COMMENT block and write to output file + comments = read_sinex_comments(sinex) + for i in comments: + out.write(f"{i}\n") + + out.write("*-------------------------------------------------------------------------------\n") # Read in the +SITE/ID block and write to file site_id = read_sinex_site_id_block(sinex) for line in site_id: out.write(f"{line}") del site_id + + out.write("*-------------------------------------------------------------------------------\n") # Read in the +SOLUTION/EPOCHS block and write to file solution_epochs = read_sinex_solution_epochs_block(sinex) @@ -893,12 +956,16 @@ def remove_matrixzeros_sinex(sinex): out.write(f"{line}") del solution_epochs + out.write("*-------------------------------------------------------------------------------\n") + # Read in the +SOLUTION/ESTIMATE block solution_estimate = read_sinex_solution_estimate_block(sinex) for line in solution_estimate: out.write(f"{line}") del solution_estimate + out.write("*-------------------------------------------------------------------------------\n") + # Read in the +SOLUTION/MATRIX_ESTIMATE block: # - Remove lines that contain only zeros solution_matrix_estimate = read_sinex_solution_matrix_estimate_block(sinex) From 78eb563cfd06c37729b08530acb5d9cf66ef9818 Mon Sep 17 00:00:00 2001 From: adeane-ga Date: Tue, 10 Jan 2023 16:03:43 +1100 Subject: [PATCH 6/6] Fix minor bug of float output in sinex header --- geodepy/gnss.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/geodepy/gnss.py b/geodepy/gnss.py index 9a94883..2b6cb8c 100644 --- a/geodepy/gnss.py +++ b/geodepy/gnss.py @@ -765,9 +765,9 @@ def remove_velocity_sinex(sinex): old_creation_time = header[15:27] creation_time = set_creation_time() header = header.replace(old_creation_time, creation_time) - old_num_params = header[60:65] - num_params = int(old_num_params) / 2 - header = header.replace(old_num_params, str(num_params)) + old_num_params = int(header[60:65]) + num_params = int(old_num_params / 2) + header = header.replace(str(old_num_params), str(num_params)) header = header.replace('V', '') out.write(header) del header