-
Notifications
You must be signed in to change notification settings - Fork 42
/
Copy pathexample_plot_md.py
105 lines (83 loc) · 2.96 KB
/
example_plot_md.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#
# example file for ppp-tools
# /~https://github.com/aewallin/ppp-tools
#
import datetime
import os
import matplotlib.pyplot as plt
import time
import numpy
import matplotlib
matplotlib.rcParams['axes.formatter.useoffset'] = False
"""
plot result of multi-day run
"""
# import ppp_rtklib
# import ppp_glab
import ppp_gpsppp
import station
import ppp_common
import jdutil
#station1 = station.mi02
station1 = station.mi06local
#station2 = station.mi05
station2 = station.mi05
#dt = datetime.datetime.utcnow()-datetime.timedelta(days=32) # 4 days ago
# first interval
#dt = datetime.datetime(2024,5,19)
#num_days = 5
# second interval
#dt = datetime.datetime(2024,5,29)
#num_days = 8
# third interval
dt = datetime.datetime(2024,5,29)
num_days = 8
current_dir = os.getcwd()
#num_days = 20
products = "final"
# mi05 - ptbb: +3.96e-16 4.13-16
# mi05 - pt10: +2.62e-16 4.095e-16
#
# mi05 - mi04: 5.43e-16
program = "gpspace"
current_dir = os.getcwd()
#day_list = []
#for n in [8, 7, 6, 5, 4, 3]:
#for n in [ 10, 9, 8, 7, 6, 5, 4, 3, 2]:
# day_list.append( datetime.datetime.utcnow()-datetime.timedelta(days=n) )
r1 = ppp_common.read_result_file(station1, dt, products, program, current_dir, num_days=num_days)
r2 = ppp_common.read_result_file(station2, dt, products, program, current_dir, num_days=num_days)
(t45, d45) = ppp_common.diff_stations( current_dir, station1, station2, dt, products, program, num_days=num_days)
mjd = jdutil.dt2mjd(t45)
d45 = numpy.array(d45)
#(t25, d25) = read_days( station.mi02, station.mi05, day_list )
#(t24, d24) = read_days( station.mi02, station.mi04, day_list )
# compute double difference
# (t_glab,d_glab) = ppp_common.diff_stations(current_dir, station1, station2, dt, products, "glab")
# (t_rtklib,d_rtklib) = ppp_common.diff_stations(current_dir, station1, station2, dt, products, "rtklib")
#(t_gpspace, d_gpspace) = ppp_common.diff_stations(current_dir, station1, station2, dt, products, "gpspace")
plt.figure()
plt.subplot(3,1,1)
plt.plot( jdutil.dt2mjd( [x.epoch for x in r1.observations]), [x.clock for x in r1.observations], label="%s" % (station1.name))
plt.ylabel('%s / ns' % station1.name)
plt.grid()
plt.subplot(3,1,2)
plt.plot( jdutil.dt2mjd( [x.epoch for x in r2.observations]), [x.clock for x in r2.observations], label="%s" % (station2.name))
plt.ylabel('%s / ns' % station2.name)
plt.grid()
plt.subplot(3,1,3)
plt.title("%s - %s receiver clock, double difference"%(station1.name, station2.name))
#plt.plot(t_glab,d_glab, label="glab")
#plt.plot(t_rtklib,d_rtklib, label="rtklib")
plt.plot(mjd, numpy.array(d45), label="%s -%s" % (station1.name, station2.name))
p = numpy.polyfit(mjd, d45, 1)
rate = p[0]
y = -1e-9*rate/(24.0*3600.0)
plt.plot(mjd, numpy.polyval(p, mjd),'r--',label='fit %.3f ns/day = %.6g'%(rate, y))
#plt.plot(t25, numpy.array(d25), label="2-5 gpspace")
#plt.plot(t24, d24, label="2-4 gpspace")
#plt.ylim((-3, 3))
plt.legend(loc="best")
plt.grid()
plt.ylabel("%s - %s / ns" %(station1.name, station2.name))
plt.show()