Index: trunk/BNC/scripts/plot_xyz.py =================================================================== --- trunk/BNC/scripts/plot_xyz.py (revision 10791) +++ trunk/BNC/scripts/plot_xyz.py (revision 10791) @@ -0,0 +1,233 @@ +#!/usr/bin/python3 + +import sys +import datetime +import matplotlib.pyplot as plt +import matplotlib.dates as mdates +import matplotlib.ticker as mticker +import numpy as np +import pathlib as pth + +import utils + +sys.dont_write_bytecode = True + +# Data Class +########################################################################################## +class Data: + def __init__(self): + self.epo0 = None + self.epo = [] + self.dNfix = [] + self.dEfix = [] + self.dUfix = [] + self.dNflt = [] + self.dEflt = [] + self.dUflt = [] + self.nFix = [] + self.conv2D = {} + + def append(self, epo, dN, dE, dU, isFix, nFix, flgReset): + self.epo.append(epo) + self.nFix.append(nFix) + if isFix: + self.dNfix.append(dN) + self.dEfix.append(dE) + self.dUfix.append(dU) + self.dNflt.append(np.nan) + self.dEflt.append(np.nan) + self.dUflt.append(np.nan) + else: + self.dNflt.append(dN) + self.dEflt.append(dE) + self.dUflt.append(dU) + self.dNfix.append(np.nan) + self.dEfix.append(np.nan) + self.dUfix.append(np.nan) + + if flgReset or self.epo0 is None: + self.epo0 = epo + + nSec = int((epo - self.epo0).total_seconds()) + if nSec not in self.conv2D: + self.conv2D[nSec] = [] + self.conv2D[nSec].append(np.sqrt(dN*dN + dE*dE)) + +# Function to generate the plot +########################################################################################## +def run_plot_client(staName, data, title, timeStr, pngFile, statFile): + cm = 1/2.54 + fig, (ax1, ax2, ax3, dummy, ax4) = plt.subplots(5, 1, figsize=(20*cm,20*cm), + gridspec_kw={'height_ratios': [1, 1, 1, 0.3, 1.5]}) + dummy.axis("off") + + ax1.set_title(title) + ax1.set_ylabel("North [m]") + ax1.tick_params(labelbottom=False) + ax1.plot(data.epo, data.dNflt, color = "crimson", linewidth = 1.0) + ax1.plot(data.epo, data.dNfix, color = "forestgreen", label = "fixed", linewidth = 1.0) + + ax2.set_ylabel("East [m]") + ax2.tick_params(labelbottom=False) + ax2.plot(data.epo, data.dEflt, color = "crimson", label = "float", linewidth = 1.0) + ax2.plot(data.epo, data.dEfix, color = "forestgreen", linewidth = 1.0) + + ax3.set_xlabel(timeStr) + ax3.set_ylabel("Height [m]") + ax3.xaxis.set_major_formatter(mdates.DateFormatter('%H:%M')) + ax3.plot(data.epo, data.dUflt, color = "crimson", linewidth = 1.0) + ax3.plot(data.epo, data.dUfix, color = "forestgreen", linewidth = 1.0) + + ax32 = ax3.twinx() + ax32.plot(data.epo, data.nFix, color = "lightskyblue", label = "% fix", linewidth = 0.4) + ax32.set_ylim(0,300) + ax32.set_yticks([0,50,100]) + ax32.legend(loc="upper right") + + for ax in [ax1, ax2, ax3]: + ax.set_ylim(-0.4, 0.4) + hh, ll = ax.get_legend_handles_labels() + if len(hh) > 0: + ax.legend(loc="upper right") + ax.grid(axis = 'y', which = 'both', color = 'grey', linestyle = '--', linewidth = 0.3) + ax.set_yticks([x / 10 for x in range(-4, 5, 2)]) + ax.yaxis.set_minor_locator(mticker.AutoMinorLocator(2)) + + maxLen = 0 + perc = {68: [], 95: []} + percMin = {10: {68 : 0.0, 95 : 0.0}, 20: {68 : 0.0, 95 : 0.0}, 30: {68 : 0.0, 95 : 0.0}} + keyMin = {10: {68 : 0.0, 95 : 0.0}, 20: {68 : 0.0, 95 : 0.0}, 30: {68 : 0.0, 95 : 0.0}} + + for key in data.conv2D.keys(): + if len(data.conv2D[key]) > maxLen: + maxLen = len(data.conv2D[key]) + + for pKey in (68, 95): + pp = np.percentile(data.conv2D[key], pKey) + perc[pKey].append(pp) + for minute in (10, 20, 30): + tMin = 60*minute + if key in range(60*minute-10, 60*minute+1): + if abs(key - tMin) < abs(keyMin[minute][pKey] - tMin): + keyMin[minute][pKey] = key + percMin[minute][pKey] = pp + + print("%8s 68: %6.3f %6.3f %6.3f 95: %6.3f %6.3f %6.3f" % + (staName, percMin[10][68], percMin[20][68], percMin[30][68], percMin[10][95], percMin[20][95], percMin[30][95])) + + for ii in range(0, maxLen): + dPos = [] + for key in data.conv2D.keys(): + if len(data.conv2D[key]) > ii: + dPos.append(data.conv2D[key][ii]) + else: + dPos.append(np.nan) + ax4.plot(data.conv2D.keys(), dPos, linewidth = 0.5) + + ax4.set_title("Convergence") + ax4.plot(data.conv2D.keys(), perc[68], linewidth = 1.0, color = "black", label = "68 %") + ax4.plot(data.conv2D.keys(), perc[95], linewidth = 2.0, color = "black", label = "95 %") + ax4.legend(loc="upper right") + ax4.set_ylim(0, 1.0) + ax4.set_xlabel("Seconds after Reset") + ax4.set_ylabel("Meters") + ax4.text( 0.1, 0.7, + "percentile after 10 min:\n68%% : %8.3f m\n95%% : %8.3f m" % (percMin[10][68], percMin[10][95]), + transform=ax4.transAxes) + ax4.text( 0.4, 0.5, + "percentile after 20 min:\n68%% : %8.3f m\n95%% : %8.3f m" % (percMin[20][68], percMin[20][95]), + transform=ax4.transAxes) + ax4.text( 0.7, 0.3, + "percentile after 30 min:\n68%% : %8.3f m\n95%% : %8.3f m" % (percMin[30][68], percMin[30][95]), + transform=ax4.transAxes) + + ### plt.show() + plt.savefig(pngFile) + plt.close() + + # Output Statistics + # ----------------- + if statFile: + with open(statFile, 'w') as outStat: + print("%s %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f" % + (fileName, + percMin[10][68], percMin[10][95], + percMin[20][68], percMin[20][95], + percMin[30][68], percMin[30][95]), file = outStat) + + +# Main Program +########################################################################################## +if __name__ == '__main__': + import sys + import argparse + + parser = argparse.ArgumentParser() + parser.add_argument("crdFile") + parser.add_argument("fileName") + parser.add_argument("--pngFile", type=str) + parser.add_argument("--statFile", type=str) + parser.add_argument("--title", type=str) + args = parser.parse_args() + + crdFile = args.crdFile + fileName = args.fileName + if args.pngFile is None: + pngFile = pth.Path(fileName).with_suffix(".png") + else : + pngFile = args.pngFile + if args.title is None: + title = fileName + else: + title = args.title + statFile = args.statFile + + dateFmt = "%Y-%m-%d_%H:%M:%S.%f" + dfCrd = utils.readcrdfile(crdFile) + station = None + + # Read Data + # --------- + data = Data() + flgReset = True + timeStr = None + with open(fileName, 'r') as inFile: + for line in inFile: + + if line.find("RESET FILTER") == 0: + flgReset = True + + elif line.find("X =") >= 0: + fields = line.split() + + dateStr = fields[0] + staName = fields[1] + xx = float(fields[4]) + yy = float(fields[9]) + zz = float(fields[14]) + isFix = (fields[32] == "fix") + if isFix: + numFix = int(fields[33]) + else: + numFix = 0 + + if station is None: + station = utils.getsta(dfCrd, staName) + if station is None: + raise Exception("Station %s not found" % staName) + ell = utils.xyz2ell((xx, yy, zz)) + + xyz = (xx - station.xx, + yy - station.yy, + zz - station.zz) + + neu = utils.xyz2neu(ell, xyz) + + time = datetime.datetime.strptime(dateStr, dateFmt) + data.append(time, neu[0], neu[1], neu[2], isFix, numFix, flgReset) + flgReset = False + + if timeStr is None: + timeStr = time.strftime("%Y-%m-%d") + + run_plot_client(staName, data, title, timeStr, pngFile, statFile) Index: trunk/BNC/scripts/utils.py =================================================================== --- trunk/BNC/scripts/utils.py (revision 10791) +++ trunk/BNC/scripts/utils.py (revision 10791) @@ -0,0 +1,270 @@ + +import sys +import pandas as pd +import numpy as np +import datetime +import colorsys +import matplotlib +from enum import Enum as Enum + +from pandas.plotting import register_matplotlib_converters +register_matplotlib_converters() + +sys.dont_write_bytecode = True + +############################################################################################## +class GenConst: + aell = 6378137.000; + fInv = 298.2572236; + rho_deg = 180.0 / np.pi; + def __init__(self): + raise TypeError("Constants class cannot be instantiated.") + + def __setattr__(self, key, value): + raise AttributeError("Constants cannot be modified.") + +############################################################################################## +class RetVal: + class Value(Enum): + SUCCESS = 0 + FAILURE = 1 + + def __init__(self): + self.value = Value.FAILURE + self.errmsg = "Uknown error" + + def is_success(self): + return self.value == VALUE.SUCCESS + + def is_failure(self): + return self.value == VALUE.FAILURE + + @staticmethod + def success(): + retval = RetVal() + retval.value = Value.SUCCESS + return retval + + @staticmethod + def failure(errmsg): + retval = RetVal() + retval.value = Value.FAILURE + retval.errmsg = errmsg + return retval + +############################################################################################## +class SatId: + def __init__(self, sys, svn): + if (sys == 'G' or + sys == 'R' or + sys == 'E' or + sys == 'C' or + sys == 'J' or + sys == 'I') : + self.sys = sys + self.svn = svn + else: + self.sys = '?' + self.svn = 0 + + def set(self, satid): + try: + if len(satid) == 0: + self.sys = '?' + self.svn = 0 + elif len(satid) == 3: + if (satid[0] == 'G' or + satid[0] == 'R' or + satid[0] == 'E' or + satid[0] == 'C' or + satid[0] == 'J' or + satid[0] == 'I') : + self.sys = satid[0] + self.svn = int(satid[1:]) + else: + RaiseException("") + else: + RaiseException("") + except: + self.sys = '?' + self.svn = 0 + + def __str__(self): + return "%s%02d" % (self.sys, self.svn) + + def __hash__(self): + return hash((self.sys, self.svn)) + + def __eq__(self,other): + return self.sys == other.sys and self.svn == other.svn + + +############################################################################################## +class StaInfo: + def __init__(self, name, xx, yy, zz): + self.name = name + self.xx = xx + self.yy = yy + self.zz = zz + + def print(self): + print("%6s %12.3f %12.3f %12.3f" % (self.name, self.xx, self.yy, self.zz)) + +############################################################################################## +def readcrdfile(filename): + df = pd.read_csv(filename, sep="\\s+", header=None, comment='#') + return df.apply(lambda row: StaInfo(row[0], row[1], row[2], row[3]), axis=1).tolist() + +############################################################################################## +def readstatfile(filename): + df = pd.read_csv(filename, sep="\\s+", header=None, comment='#') + retval= {} + for row in df.itertuples(): + ##retval[row._1[0:4]] = [row._2, row._3, row._4, row._5, row._6, row._7] + retval[row._1] = [row._2, row._3, row._4, row._5, row._6, row._7] + return retval + +############################################################################################## +def getsta(stalist, name): + for sta in stalist: + if sta.name == name: + return sta + return None + +############################################################################################## +def xyz2ell(xyz): + bell = GenConst.aell*(1.0-1.0/GenConst.fInv) + e2 = (GenConst.aell*GenConst.aell-bell*bell)/(GenConst.aell*GenConst.aell) + e2c = (GenConst.aell*GenConst.aell-bell*bell)/(bell*bell) + + if (xyz[0] == 0.0 and xyz[1] == 0.0 and xyz[2] == 0.0): + return [] + + ss = np.sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]) ; + zps = xyz[2]/ss ; + theta = np.arctan( (xyz[2]*GenConst.aell) / (ss*bell) ); + sin3 = np.sin(theta) * np.sin(theta) * np.sin(theta); + cos3 = np.cos(theta) * np.cos(theta) * np.cos(theta); + + ell = [0, 0, 0]; + ell[0] = np.arctan( (xyz[2] + e2c * bell * sin3) / (ss - e2 * GenConst.aell * cos3) ); + ell[1] = np.arctan2(xyz[1],xyz[0]) ; + nn = GenConst.aell/np.sqrt(1.0-e2*np.sin(ell[0])*np.sin(ell[0])) ; + ell[2] = ss / np.cos(ell[0]) - nn; + + for ii in range(10): + nn = GenConst.aell/np.sqrt(1.0-e2*np.sin(ell[0])*np.sin(ell[0])) + hOld = ell[2] + phiOld = ell[0] + ell[2] = ss/np.cos(ell[0])-nn + ell[0] = np.arctan(zps/(1.0-e2*nn/(nn+ell[2]))) + if ( abs(phiOld-ell[0]) <= 1.0e-11 and abs(hOld-ell[2]) <= 1.0e-5 ): + return ell; + +############################################################################################## +def ell2xyz(ell): + bell = GenConst.aell*(1.0-1.0/GenConst.fInv) + + e2 = (GenConst.aell*GenConst.aell-bell*bell)/(GenConst.aell*GenConst.aell) + sinPhi = np.sin(ell[0]) + cosPhi = np.sqrt(1.0-sinPhi*sinPhi) + nn = GenConst.aell/np.sqrt(1.0-e2*sinPhi*sinPhi) + + return [ (nn + ell[2]) * cosPhi * np.cos(ell[1]), + (nn + ell[2]) * cosPhi * np.sin(ell[1]), + (nn*(1.0-e2)+ ell[2]) * sinPhi ] + +############################################################################################## +def mjd2datetime(mjd): + EPOCH0 = dtatetime.datetime(1980, 1, 6) + MJD0 = 44244.0 + return EPOCH0 + timedelta(days=mjd-MJD0) + +############################################################################################## +def gpswkdow2datetime(gpswk, dow): + EPOCH0 = datetime.datetime(1980, 1, 6) + return EPOCH0 + timedelta(days=gpswk*7+dow) + +############################################################################################## +def jacobiXyz2neu(ell): + sinPhi = np.sin(ell[0]) + cosPhi = np.cos(ell[0]) + sinLam = np.sin(ell[1]) + cosLam = np.cos(ell[1]) + + return np.array([[-sinPhi*cosLam, -sinPhi*sinLam, cosPhi], + [-sinLam, cosLam, 0.0 ], + [cosPhi*cosLam, cosPhi*sinLam, sinPhi]]) + +############################################################################################## +def xyz2neu(ell, xyz): + return jacobiXyz2neu(ell).dot(xyz) + +########################################################################################## +def modColor(cc, light, satur): + hh, ll, ss = colorsys.rgb_to_hls(*matplotlib.colors.to_rgb(cc)) + return colorsys.hls_to_rgb(hh, light*ll, satur) + +############################################################################################## +def currentgpswk(): + now = datetime.now() + return gpswk(now) + +############################################################################################## +def currentgpswkdow(): + now = datetime.now() + return gpswkdow(now) + +############################################################################################## +def getcurrentsp3file(localoutdir): + tt = datetime.now(tz=timezone.utc) - timedelta(hours=2) + tnow = datetime.now(tz=timezone.utc) + + while True: + gpsw, dow = gpswkdow(tt) + hours = tt.hour + minutes = int(tt.minute / 15) * 15 + sp3filenameLocal = "%s/vmcpp%04d%d_%02d%02d.sp3" % (localoutdir, gpsw, dow, hours, minutes) + sp3filenameSftp = "%04d/vmcpp%04d%d_%02d%02d.sp3" % (gpsw, gpsw, dow, hours, minutes) + + if not os.path.exists(sp3filenameLocal): + command = "echo -e \"lcd %s\\\\n get %s\" | sftp -i ~/.ssh/ocds_archive.rsa ocdsarchive@napeos.correction-products.veripos.com" % (localoutdir, sp3filenameSftp) + print("Try to download file %s, cmd %s" % (sp3filenameSftp, command)) + return_code = subprocess.run(command, shell=True, stdout = subprocess.DEVNULL, stderr = subprocess.DEVNULL) + ###return_code = subprocess.run(command, shell=True) + + if os.path.exists(sp3filenameLocal): + print(" download succeed") + else: + print(" download failed") + else: + print("File %s exists already\n" % (sp3filenameLocal)) + + tt = tt + timedelta(minutes=15) + if tt > tnow: + break + +############################################################################################## +def deleteoldfiles(extension: str, nhours: float, dirname: str): + if not extension.startswith('.'): + extension = '.' + extension + + dirpath = Path(dirname) + if not dirpath.is_dir(): + print(f"Error: Directory not found: {dirname}") + return + + # Compute threshold time in seconds + now = time.time() + age_limit = now - nhours * 3600 + + deleted_count = 0 + for file in dirpath.glob(f"*{extension}"): + try: + mtime = file.stat().st_mtime + if mtime < age_limit: + file.unlink() + deleted_count += 1 + print(f"Deleted: {file}") + except Exception as e: + print(f"Failed to delete {file}: {e}")