From f651b81d6aaf92b086e28d45e30f38b6ac9225fd Mon Sep 17 00:00:00 2001 From: Craig Harrison Date: Sun, 17 Jul 2022 08:55:37 +1000 Subject: [PATCH] Upper triangular matrices now supported in remove_stns_sinex() --- .gitignore | 11 ++++++++++ geodepy/gnss.py | 54 +++++++++++++++++++++++++++++-------------------- 2 files changed, 43 insertions(+), 22 deletions(-) diff --git a/.gitignore b/.gitignore index be736c0..4e7e7d4 100644 --- a/.gitignore +++ b/.gitignore @@ -34,3 +34,14 @@ env /geodepy/XVSOLFIN_20220326.SNX /geodepy/exciseStationsAPREF.py /geodepy/gda2020.dat +/geodepy/something_list_snx_blks.py +/geodepy/try_remove_stns.py +/geodepy/17143.U.SNX.AUS +/geodepy/17143.L.SNX.AUS +/geodepy/try_remove_stns.py +/geodepy/SNXEPO.CON.SNX +/geodepy/SNXEPO.NONCON.SNX +/geodepy/SNXEPO.SNX.CON.AUS +/geodepy/SNXEPO.SNX.NONCON.AUS +/geodepy/output.l.snx +/geodepy/output.u.snx diff --git a/geodepy/gnss.py b/geodepy/gnss.py index 001e3b6..e5b510e 100644 --- a/geodepy/gnss.py +++ b/geodepy/gnss.py @@ -7,7 +7,6 @@ In Development """ -import sys from datetime import datetime from numpy import zeros from geodepy.angles import DMSAngle @@ -44,6 +43,7 @@ def print_sinex_comments(file): if line.startswith('-FILE/COMMENT'): go = False + def set_creation_time(): """This function sets the creation time, in format YY:DDD:SSSSS, for use in the SINEX header line @@ -451,7 +451,7 @@ def read_sinex_header_block(sinex): block = [] with open(sinex, 'r') as f: - line = f.readline() + next(f) line = f.readline() while line: block.append(line) @@ -561,9 +561,6 @@ def remove_stns_sinex(sinex, sites): :return: SINEX file output.snx """ - # The block separator - separator = '*' + '-' * 79 + '\n' - # Open the output file with open('output.snx', 'w') as out: @@ -573,7 +570,7 @@ def remove_stns_sinex(sinex, sites): old_creation_time = header[15:27] creation_time = set_creation_time() header = header.replace(old_creation_time, creation_time) - old_num_params = int(header[60:65]) + old_num_params = header[60:65] if header[70:71] == 'V': num_stn_params = 6 else: @@ -585,11 +582,10 @@ def remove_stns_sinex(sinex, sites): if site in sites: num_stns_to_remove += 1 del solution_epochs - num_params = old_num_params - num_stn_params * num_stns_to_remove + num_params = int(old_num_params) - num_stn_params * num_stns_to_remove num_params = '{:05d}'.format(num_params) header = header.replace(str(old_num_params), str(num_params)) out.write(header) - out.write(separator) # Read in the site ID block and write out the sites not being removed site_id = read_sinex_site_id_block(sinex) @@ -602,7 +598,6 @@ def remove_stns_sinex(sinex, sites): if site not in sites: out.write(line) del site_id - out.write(separator) # Read in the solution epochs block and write out the epochs of the # sites not being removed @@ -616,7 +611,6 @@ def remove_stns_sinex(sinex, sites): if site not in sites: out.write(line) del solution_epochs - out.write(separator) # Read in the solution estimate block and write out the estimates of # the sites not being removed @@ -638,15 +632,19 @@ def remove_stns_sinex(sinex, sites): line = ' ' + number + line[6:] out.write(line) del solution_estimate - out.write(separator) # Read in the matrix estimate block and write out minus the sites # being removed vcv = {} 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' out.write(solution_matrix_estimate[0]) - out.write(solution_matrix_estimate[1]) + if solution_matrix_estimate[1].startswith('*'): + out.write(solution_matrix_estimate[1]) for line in solution_matrix_estimate: if line.startswith(' '): cols = line.split() @@ -664,24 +662,36 @@ def remove_stns_sinex(sinex, sites): for i in range(1, len(vcv)+1): if i not in skip: sub_row += 1 - for j in range(i): - if j+1 not in skip: - try: - sub_vcv[str(sub_row)].append(vcv[str(i)][j]) - except KeyError: - sub_vcv[str(sub_row)] = [] - sub_vcv[str(sub_row)].append(vcv[str(i)][j]) + if matrix == 'lower': + for j in range(i): + if j+1 not in skip: + try: + sub_vcv[str(sub_row)].append(vcv[str(i)][j]) + except KeyError: + sub_vcv[str(sub_row)] = [] + sub_vcv[str(sub_row)].append(vcv[str(i)][j]) + if matrix == 'upper': + for j in range(len(vcv)-(i-1)): + if j+i not in skip: + try: + sub_vcv[str(sub_row)].append(vcv[str(i)][j]) + except KeyError: + sub_vcv[str(sub_row)] = [] + sub_vcv[str(sub_row)].append(vcv[str(i)][j]) for i in range(1, len(sub_vcv)+1): para1 = '{:5d}'.format(i) - j = -2 + if matrix == 'lower': + j = -2 + elif matrix == 'upper': + j = i - 3 while sub_vcv[str(i)]: j += 3 para2 = '{:5d}'.format(j) line = ' ' + para1 + ' ' + para2 n = min([3, len(sub_vcv[str(i)])]) for k in range(n): - val = sub_vcv[str(i)].pop(0) - line += ' ' + val + val = '{:21.14e}'.format(float(sub_vcv[str(i)].pop(0))) + line += ' ' + str(val) out.write(line + '\n') out.write(block_end)