source: ntrip/trunk/BNC/scripts/utils.py

Last change on this file was 10791, checked in by mervart, 4 weeks ago

BNC Multifrequency and PPPAR Client (initial version)
File size: 18.7 KB
Line 
1
2import sys
3import pandas as pd
4import numpy as np
5import datetime
6import colorsys
7import matplotlib
8from enum import Enum as Enum
9
10from pandas.plotting import register_matplotlib_converters
11register_matplotlib_converters()
12
13sys.dont_write_bytecode = True
14
15##############################################################################################
16class GenConst:
17 aell = 6378137.000;
18 fInv = 298.2572236;
19 rho_deg = 180.0 / np.pi;
20 def __init__(self):
21 raise TypeError("Constants class cannot be instantiated.")
22
23 def __setattr__(self, key, value):
24 raise AttributeError("Constants cannot be modified.")
25
26##############################################################################################
27class RetVal:
28 class Value(Enum):
29 SUCCESS = 0
30 FAILURE = 1
31
32 def __init__(self):
33 self.value = Value.FAILURE
34 self.errmsg = "Uknown error"
35
36 def is_success(self):
37 return self.value == VALUE.SUCCESS
38
39 def is_failure(self):
40 return self.value == VALUE.FAILURE
41
42 @staticmethod
43 def success():
44 retval = RetVal()
45 retval.value = Value.SUCCESS
46 return retval
47
48 @staticmethod
49 def failure(errmsg):
50 retval = RetVal()
51 retval.value = Value.FAILURE
52 retval.errmsg = errmsg
53 return retval
54
55##############################################################################################
56class SatId:
57 def __init__(self, sys, svn):
58 if (sys == 'G' or
59 sys == 'R' or
60 sys == 'E' or
61 sys == 'C' or
62 sys == 'J' or
63 sys == 'I') :
64 self.sys = sys
65 self.svn = svn
66 else:
67 self.sys = '?'
68 self.svn = 0
69
70 def set(self, satid):
71 try:
72 if len(satid) == 0:
73 self.sys = '?'
74 self.svn = 0
75 elif len(satid) == 3:
76 if (satid[0] == 'G' or
77 satid[0] == 'R' or
78 satid[0] == 'E' or
79 satid[0] == 'C' or
80 satid[0] == 'J' or
81 satid[0] == 'I') :
82 self.sys = satid[0]
83 self.svn = int(satid[1:])
84 else:
85 RaiseException("")
86 else:
87 RaiseException("")
88 except:
89 self.sys = '?'
90 self.svn = 0
91
92 def __str__(self):
93 return "%s%02d" % (self.sys, self.svn)
94
95 def __hash__(self):
96 return hash((self.sys, self.svn))
97
98 def __eq__(self,other):
99 return self.sys == other.sys and self.svn == other.svn
100
101
102##############################################################################################
103class StaInfo:
104 def __init__(self, name, xx, yy, zz):
105 self.name = name
106 self.xx = xx
107 self.yy = yy
108 self.zz = zz
109
110 def print(self):
111 print("%6s %12.3f %12.3f %12.3f" % (self.name, self.xx, self.yy, self.zz))
112
113##############################################################################################
114def readcrdfile(filename):
115 df = pd.read_csv(filename, sep="\\s+", header=None, comment='#')
116 return df.apply(lambda row: StaInfo(row[0], row[1], row[2], row[3]), axis=1).tolist()
117
118##############################################################################################
119def readstatfile(filename):
120 df = pd.read_csv(filename, sep="\\s+", header=None, comment='#')
121 retval= {}
122 for row in df.itertuples():
123 ##retval[row._1[0:4]] = [row._2, row._3, row._4, row._5, row._6, row._7]
124 retval[row._1] = [row._2, row._3, row._4, row._5, row._6, row._7]
125 return retval
126
127##############################################################################################
128def getsta(stalist, name):
129 for sta in stalist:
130 if sta.name == name:
131 return sta
132 return None
133
134##############################################################################################
135def xyz2ell(xyz):
136 bell = GenConst.aell*(1.0-1.0/GenConst.fInv)
137 e2 = (GenConst.aell*GenConst.aell-bell*bell)/(GenConst.aell*GenConst.aell)
138 e2c = (GenConst.aell*GenConst.aell-bell*bell)/(bell*bell)
139
140 if (xyz[0] == 0.0 and xyz[1] == 0.0 and xyz[2] == 0.0):
141 return []
142
143 ss = np.sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]) ;
144 zps = xyz[2]/ss ;
145 theta = np.arctan( (xyz[2]*GenConst.aell) / (ss*bell) );
146 sin3 = np.sin(theta) * np.sin(theta) * np.sin(theta);
147 cos3 = np.cos(theta) * np.cos(theta) * np.cos(theta);
148
149 ell = [0, 0, 0];
150 ell[0] = np.arctan( (xyz[2] + e2c * bell * sin3) / (ss - e2 * GenConst.aell * cos3) );
151 ell[1] = np.arctan2(xyz[1],xyz[0]) ;
152 nn = GenConst.aell/np.sqrt(1.0-e2*np.sin(ell[0])*np.sin(ell[0])) ;
153 ell[2] = ss / np.cos(ell[0]) - nn;
154
155 for ii in range(10):
156 nn = GenConst.aell/np.sqrt(1.0-e2*np.sin(ell[0])*np.sin(ell[0]))
157 hOld = ell[2]
158 phiOld = ell[0]
159 ell[2] = ss/np.cos(ell[0])-nn
160 ell[0] = np.arctan(zps/(1.0-e2*nn/(nn+ell[2])))
161 if ( abs(phiOld-ell[0]) <= 1.0e-11 and abs(hOld-ell[2]) <= 1.0e-5 ):
162 return ell;
163
164##############################################################################################
165def ell2xyz(ell):
166 bell = GenConst.aell*(1.0-1.0/GenConst.fInv)
167
168 e2 = (GenConst.aell*GenConst.aell-bell*bell)/(GenConst.aell*GenConst.aell)
169 sinPhi = np.sin(ell[0])
170 cosPhi = np.sqrt(1.0-sinPhi*sinPhi)
171 nn = GenConst.aell/np.sqrt(1.0-e2*sinPhi*sinPhi)
172
173 return [ (nn + ell[2]) * cosPhi * np.cos(ell[1]),
174 (nn + ell[2]) * cosPhi * np.sin(ell[1]),
175 (nn*(1.0-e2)+ ell[2]) * sinPhi ]
176
177##############################################################################################
178def mjd2datetime(mjd):
179 EPOCH0 = dtatetime.datetime(1980, 1, 6)
180 MJD0 = 44244.0
181 return EPOCH0 + timedelta(days=mjd-MJD0)
182
183##############################################################################################
184def gpswkdow2datetime(gpswk, dow):
185 EPOCH0 = datetime.datetime(1980, 1, 6)
186 return EPOCH0 + timedelta(days=gpswk*7+dow)
187
188##############################################################################################
189def jacobiXyz2neu(ell):
190 sinPhi = np.sin(ell[0])
191 cosPhi = np.cos(ell[0])
192 sinLam = np.sin(ell[1])
193 cosLam = np.cos(ell[1])
194
195 return np.array([[-sinPhi*cosLam, -sinPhi*sinLam, cosPhi],
196 [-sinLam, cosLam, 0.0 ],
197 [cosPhi*cosLam, cosPhi*sinLam, sinPhi]])
198
199##############################################################################################
200def xyz2neu(ell, xyz):
201 return jacobiXyz2neu(ell).dot(xyz)
202
203##########################################################################################
204def modColor(cc, light, satur):
205 hh, ll, ss = colorsys.rgb_to_hls(*matplotlib.colors.to_rgb(cc))
206 return colorsys.hls_to_rgb(hh, light*ll, satur)
207
208##############################################################################################
209def currentgpswk():
210 now = datetime.now()
211 return gpswk(now)
212
213##############################################################################################
214def currentgpswkdow():
215 now = datetime.now()
216 return gpswkdow(now)
217
218##############################################################################################
219def getcurrentsp3file(localoutdir):
220 tt = datetime.now(tz=timezone.utc) - timedelta(hours=2)
221 tnow = datetime.now(tz=timezone.utc)
222
223 while True:
224 gpsw, dow = gpswkdow(tt)
225 hours = tt.hour
226 minutes = int(tt.minute / 15) * 15
227 sp3filenameLocal = "%s/vmcpp%04d%d_%02d%02d.sp3" % (localoutdir, gpsw, dow, hours, minutes)
228 sp3filenameSftp = "%04d/vmcpp%04d%d_%02d%02d.sp3" % (gpsw, gpsw, dow, hours, minutes)
229
230 if not os.path.exists(sp3filenameLocal):
231 command = "echo -e \"lcd %s\\\\n get %s\" | sftp -i ~/.ssh/ocds_archive.rsa ocdsarchive@napeos.correction-products.veripos.com" % (localoutdir, sp3filenameSftp)
232 print("Try to download file %s, cmd %s" % (sp3filenameSftp, command))
233 return_code = subprocess.run(command, shell=True, stdout = subprocess.DEVNULL, stderr = subprocess.DEVNULL)
234 ###return_code = subprocess.run(command, shell=True)
235
236 if os.path.exists(sp3filenameLocal):
237 print(" download succeed")
238 else:
239 print(" download failed")
240 else:
241 print("File %s exists already\n" % (sp3filenameLocal))
242
243 tt = tt + timedelta(minutes=15)
244 if tt > tnow:
245 break
246
247##############################################################################################
248def deleteoldfiles(extension: str, nhours: float, dirname: str):
249 if not extension.startswith('.'):
250 extension = '.' + extension
251
252 dirpath = Path(dirname)
253 if not dirpath.is_dir():
254 print(f"Error: Directory not found: {dirname}")
255 return
256
257 # Compute threshold time in seconds
258 now = time.time()
259 age_limit = now - nhours * 3600
260
261 deleted_count = 0
262 for file in dirpath.glob(f"*{extension}"):
263 try:
264 mtime = file.stat().st_mtime
265 if mtime < age_limit:
266 file.unlink()
267 deleted_count += 1
268 print(f"Deleted: {file}")
269 except Exception as e:
270 print(f"Failed to delete {file}: {e}")
Note: See TracBrowser for help on using the repository browser.