Skip to content

Commit a654acf

Browse files
harry093harry093
authored
Merge pull request #149 from GeoscienceAustralia/AD/nadj_workflow
Ad/nadj workflow
2 parents 583f8fb + 78eb563 commit a654acf

File tree

1 file changed

+290
-1
lines changed

1 file changed

+290
-1
lines changed

‎geodepy/gnss.py‎

Lines changed: 290 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,8 +8,9 @@
88
"""
99

1010
from datetime import datetime
11-
from numpy import zeros
11+
from numpy import zeros, delete
1212
from geodepy.angles import DMSAngle
13+
import re
1314

1415

1516
def list_sinex_blocks(file):
@@ -43,6 +44,28 @@ def print_sinex_comments(file):
4344
if line.startswith('-FILE/COMMENT'):
4445
go = False
4546

47+
def read_sinex_comments(file):
48+
"""This function reads comments in a SINEX file.
49+
50+
:param str file: the input SINEX file
51+
:return: comments block
52+
:rtype: list of strings
53+
"""
54+
comments = []
55+
go = False
56+
with open(file) as f:
57+
for line in f.readlines():
58+
if line.startswith('+FILE/COMMENT'):
59+
go = True
60+
if go:
61+
comments.append(line.strip())
62+
if line.startswith('-FILE/COMMENT'):
63+
go = False
64+
65+
comments.insert(-1,f"* File created by Geodepy.gnss.py at {datetime.now().strftime('%d-%m-%Y, %H:%M')}")
66+
67+
return comments
68+
4669

4770
def set_creation_time():
4871
"""This function sets the creation time, in format YY:DDD:SSSSS, for use
@@ -587,6 +610,15 @@ def remove_stns_sinex(sinex, sites):
587610
header = header.replace(str(old_num_params), str(num_params))
588611
out.write(header)
589612

613+
out.write("*-------------------------------------------------------------------------------\n")
614+
615+
# Read in the +FILE/COMMENT block and write to output file
616+
comments = read_sinex_comments(sinex)
617+
for i in comments:
618+
out.write(f"{i}\n")
619+
620+
out.write("*-------------------------------------------------------------------------------\n")
621+
590622
# Read in the site ID block and write out the sites not being removed
591623
site_id = read_sinex_site_id_block(sinex)
592624
for line in site_id:
@@ -599,6 +631,8 @@ def remove_stns_sinex(sinex, sites):
599631
out.write(line)
600632
del site_id
601633

634+
out.write("*-------------------------------------------------------------------------------\n")
635+
602636
# Read in the solution epochs block and write out the epochs of the
603637
# sites not being removed
604638
solution_epochs = read_sinex_solution_epochs_block(sinex)
@@ -612,6 +646,8 @@ def remove_stns_sinex(sinex, sites):
612646
out.write(line)
613647
del solution_epochs
614648

649+
out.write("*-------------------------------------------------------------------------------\n")
650+
615651
# Read in the solution estimate block and write out the estimates of
616652
# the sites not being removed
617653
skip = []
@@ -633,6 +669,8 @@ def remove_stns_sinex(sinex, sites):
633669
out.write(line)
634670
del solution_estimate
635671

672+
out.write("*-------------------------------------------------------------------------------\n")
673+
636674
# Read in the matrix estimate block and write out minus the sites
637675
# being removed
638676
vcv = {}
@@ -699,3 +737,254 @@ def remove_stns_sinex(sinex, sites):
699737
out.write('%ENDSNX\n')
700738

701739
return
740+
741+
def remove_velocity_sinex(sinex):
742+
"""This function reads in a SINEX file and removes the
743+
velocity parameters, including the zeros of the SOLUTION/MATRIX_ESTIMATE block,
744+
745+
:param str sinex: input SINEX file
746+
return: SINEX file output.snx
747+
"""
748+
749+
# From header, confirm that the SINEX has velocity parameters
750+
header = read_sinex_header_line(sinex)
751+
if header[70:71] == 'V':
752+
pass
753+
else:
754+
print("Not removing velocities because SINEX file does not have velocity parameters.")
755+
exit()
756+
757+
# Open the output file
758+
with open('output.snx', 'w') as out:
759+
760+
# With header line:
761+
# - update the creation time
762+
# - update number of parameter estimates
763+
# - remove 'V' from parameter list
764+
# - then write to file
765+
old_creation_time = header[15:27]
766+
creation_time = set_creation_time()
767+
header = header.replace(old_creation_time, creation_time)
768+
old_num_params = int(header[60:65])
769+
num_params = int(old_num_params / 2)
770+
header = header.replace(str(old_num_params), str(num_params))
771+
header = header.replace('V', '')
772+
out.write(header)
773+
del header
774+
775+
out.write("*-------------------------------------------------------------------------------\n")
776+
777+
# Read in the +FILE/COMMENT block and write to output file
778+
comments = read_sinex_comments(sinex)
779+
for i in comments:
780+
out.write(f"{i}\n")
781+
782+
out.write("*-------------------------------------------------------------------------------\n")
783+
784+
# Read in the +SITE/ID block and write to file
785+
site_id = read_sinex_site_id_block(sinex)
786+
for line in site_id:
787+
out.write(f"{line}")
788+
del site_id
789+
790+
out.write("*-------------------------------------------------------------------------------\n")
791+
792+
# Read in the +SOLUTION/EPOCHS block and write to file
793+
# - also collect count on number of sites for use later
794+
numSites = 0
795+
solution_epochs = read_sinex_solution_epochs_block(sinex)
796+
for line in solution_epochs:
797+
out.write(f"{line}")
798+
if line[0]!="+" and line[0]!="*" and line[0]!="-":
799+
numSites+=1
800+
del solution_epochs
801+
802+
out.write("*-------------------------------------------------------------------------------\n")
803+
804+
# Read in the +SOLUTION/ESTIMATE block:
805+
# - gather velocity indices
806+
# - remove velocity rows
807+
# - change indices to account for removed velocities
808+
# - write to file
809+
vel_indices = []
810+
estimate_number = 0
811+
solution_estimate = read_sinex_solution_estimate_block(sinex)
812+
for line in solution_estimate:
813+
if line[7:10]=="VEL":
814+
vel_indices.append(int(line[0:6]))
815+
elif line[0]=="+" or line[0]=="*" or line[0]=="-":
816+
out.write(f"{line}")
817+
else:
818+
estimate_number+=1
819+
number = '{:5d}'.format(estimate_number)
820+
line = ' ' + number + line[6:]
821+
out.write(f"{line}")
822+
del solution_estimate
823+
824+
out.write("*-------------------------------------------------------------------------------\n")
825+
826+
# Read in the +SOLUTION/MATRIX_ESTIMATE block:
827+
# - identify if matrix is upper or lower triangle form
828+
# - form full matrix
829+
# - remove all rows/columns associated with velocity
830+
# - write to file with new matrix indices (para1, para2)
831+
solution_matrix_estimate = read_sinex_solution_matrix_estimate_block(sinex)
832+
if solution_matrix_estimate[0][26:27] == 'L':
833+
matrix = 'lower'
834+
elif solution_matrix_estimate[0][26:27] == 'U':
835+
matrix = 'upper'
836+
# Size of original matrix
837+
matrix_dim = numSites*6
838+
# Setup matrix of zeros
839+
Q = zeros((matrix_dim, matrix_dim))
840+
# Write initial comment line(s), save last comment line, and form matrix
841+
for line in solution_matrix_estimate:
842+
if line[0]=="+":
843+
out.write(f"{line}")
844+
continue
845+
if line[0]=="*":
846+
out.write(f"{line}")
847+
continue
848+
if line[0]=="-":
849+
block_end = line
850+
else:
851+
lineMAT = re.split('\s+', line)
852+
lineMAT = list(filter(None, lineMAT))
853+
row = int(lineMAT[0])
854+
col = int(lineMAT[1])
855+
if len(lineMAT) >= 3:
856+
q1 = float(lineMAT[2])
857+
Q[row-1, col-1] = q1
858+
Q[col-1, row-1] = q1
859+
if len(lineMAT) >= 4:
860+
q2 = float(lineMAT[3])
861+
Q[row-1, col] = q2
862+
Q[col, row-1] = q2
863+
if len(lineMAT) >= 5:
864+
q3 = float(lineMAT[4])
865+
Q[row-1, col+1] = q3
866+
Q[col+1, row-1] = q3
867+
# Remove velocity rows/columns from matrix
868+
num_removed = 0
869+
for i in vel_indices:
870+
Q = delete(Q, i-1-num_removed, 0)
871+
Q = delete(Q, i-1-num_removed, 1)
872+
num_removed += 1
873+
# Write matrix to SINEX (lower triangle)
874+
if matrix == 'lower':
875+
i = 0
876+
j = 0
877+
while i < len(Q):
878+
while j <= i:
879+
out.write(f" {i+1:5d} {j+1:5d} {Q[i,j]:21.14E} ")
880+
j += 1
881+
if j <= i:
882+
out.write(f"{Q[i,j]:21.14E} ")
883+
j += 1
884+
if j <= i:
885+
out.write(f"{Q[i,j]:21.14E}")
886+
j += 1
887+
out.write(" \n")
888+
j = 0
889+
i += 1
890+
# Write matrix to SINEX (upper triangle)
891+
if matrix == 'upper':
892+
j = 0
893+
for i in range(len(Q)):
894+
j = i
895+
while j < len(Q):
896+
out.write(f" {i+1:5d} {j+1:5d} {Q[i,j]:21.14E} ")
897+
j += 1
898+
if j < len(Q):
899+
out.write(f"{Q[i,j]:21.14E} ")
900+
j += 1
901+
if j < len(Q):
902+
out.write(f"{Q[i,j]:21.14E}")
903+
j += 1
904+
out.write(" \n")
905+
# Write out end of block line, and delete large variables
906+
out.write(block_end)
907+
del solution_matrix_estimate
908+
del Q
909+
910+
# Write out the trailer line
911+
out.write('%ENDSNX\n')
912+
913+
return
914+
915+
def remove_matrixzeros_sinex(sinex):
916+
"""This function reads in a SINEX file and removes the
917+
zeros from the SOLUTION/MATRIX_ESTIMATE block only.
918+
919+
:param str sinex: input SINEX file
920+
return: SINEX file output.snx
921+
"""
922+
923+
# Open the output file
924+
with open('output.snx', 'w') as out:
925+
926+
# With header line:
927+
# - update the creation time
928+
# - then write to file
929+
header = read_sinex_header_line(sinex)
930+
old_creation_time = header[15:27]
931+
creation_time = set_creation_time()
932+
header = header.replace(old_creation_time, creation_time)
933+
out.write(header)
934+
del header
935+
936+
out.write("*-------------------------------------------------------------------------------\n")
937+
938+
# Read in the +FILE/COMMENT block and write to output file
939+
comments = read_sinex_comments(sinex)
940+
for i in comments:
941+
out.write(f"{i}\n")
942+
943+
out.write("*-------------------------------------------------------------------------------\n")
944+
945+
# Read in the +SITE/ID block and write to file
946+
site_id = read_sinex_site_id_block(sinex)
947+
for line in site_id:
948+
out.write(f"{line}")
949+
del site_id
950+
951+
out.write("*-------------------------------------------------------------------------------\n")
952+
953+
# Read in the +SOLUTION/EPOCHS block and write to file
954+
solution_epochs = read_sinex_solution_epochs_block(sinex)
955+
for line in solution_epochs:
956+
out.write(f"{line}")
957+
del solution_epochs
958+
959+
out.write("*-------------------------------------------------------------------------------\n")
960+
961+
# Read in the +SOLUTION/ESTIMATE block
962+
solution_estimate = read_sinex_solution_estimate_block(sinex)
963+
for line in solution_estimate:
964+
out.write(f"{line}")
965+
del solution_estimate
966+
967+
out.write("*-------------------------------------------------------------------------------\n")
968+
969+
# Read in the +SOLUTION/MATRIX_ESTIMATE block:
970+
# - Remove lines that contain only zeros
971+
solution_matrix_estimate = read_sinex_solution_matrix_estimate_block(sinex)
972+
for line in solution_matrix_estimate:
973+
col = line.split()
974+
numCol = len(col)
975+
if numCol==3:
976+
if col[2]=="0.00000000000000e+00":
977+
continue
978+
if numCol==4:
979+
if col[2]=="0.00000000000000e+00" and col[3]=="0.00000000000000e+00":
980+
continue
981+
if numCol==5:
982+
if col[2]=="0.00000000000000e+00" and col[3]=="0.00000000000000e+00" and col[4]=="0.00000000000000e+00":
983+
continue
984+
out.write(line)
985+
del solution_matrix_estimate
986+
987+
# Write out the trailer line
988+
out.write('%ENDSNX\n')
989+
990+
return

0 commit comments

Comments
 (0)