Skip to content

Commit

Permalink
Merge pull request #149 from GeoscienceAustralia/AD/nadj_workflow
Browse files Browse the repository at this point in the history 
Ad/nadj workflow
  • Loading branch information
@harry093
harry093 committed Apr 14, 2023
2 parents 583f8fb + 78eb563 commit a654acf
Showing 1 changed file with 290 additions and 1 deletion.
291 changes: 290 additions & 1 deletion geodepy/gnss.py
View file
Expand Up @@ -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):
Expand Down Expand Up @@ -43,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
Expand Down Expand Up @@ -587,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:
Expand All @@ -599,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)
Expand All @@ -612,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 = []
Expand All @@ -633,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 = {}
Expand Down Expand Up @@ -699,3 +737,254 @@ 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 SOLUTION/MATRIX_ESTIMATE block,
: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
# - remove 'V' from parameter list
# - 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 = 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

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
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

out.write("*-------------------------------------------------------------------------------\n")

# 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

out.write("*-------------------------------------------------------------------------------\n")

# Read in the +SOLUTION/MATRIX_ESTIMATE block:
# - identify if matrix is upper or lower triangle form
# - 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)
if solution_matrix_estimate[0][26:27] == 'L':
matrix = 'lower'
elif solution_matrix_estimate[0][26:27] == 'U':
matrix = 'upper'
# 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 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 = line
else:
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(lineMAT) >= 4:
q2 = float(lineMAT[3])
Q[row-1, col] = q2
Q[col, row-1] = q2
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 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 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 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)
del solution_matrix_estimate
del Q

# Write out the trailer line
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
"""

# Open the output file
with open('output.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

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)
for line in solution_epochs:
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)
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

0 comments on commit a654acf

Please sign in to comment.