from datetime import datetime, timedelta import numpy as np from pathlib import Path working_dir = Path(__file__).absolute().parent if '__file__' in dir() else Path.cwd() def to_datetime(dtlist): return datetime(int(dtlist[0]), int(dtlist[1]), int(dtlist[2]),int(dtlist[3]),int(dtlist[4]))+timedelta(seconds=float(dtlist[5])) def dt2float(dt,dt0): return ((dt - dt0) / timedelta(seconds=1)).astype('float') def get_sat_freq(sat): GPRN = int(sat[1:]) if sat[0] == 'G': #GPS FL1 = 1.57542E9 FL2 = 1.22760E9 elif sat[0] == 'R': #GLONASS Kfrq = [1, -4, 5, 6, 1, -4, 5, 6, -2, -7, 0, -1, -2, -7, 0, -1, 4, -3, 3, 2, 4, -3, 3, 2] baseFL1 = 1.602E9 baseFL2 = 1.246E9 DFRQ1_GLO = 5.625E5 DFRQ2_GLO = 4.375E5 FL1 = baseFL1 + Kfrq[GPRN - 1] * DFRQ1_GLO FL2 = baseFL2 + Kfrq[GPRN - 1] * DFRQ2_GLO GPRN = GPRN + 100 elif sat[0] == 'E': #GALILEO FL1 = 1.57542E9 FL2 = 1.20714E9 GPRN = GPRN + 200 elif sat[0] == 'C': # BEIDOU FL1 = 1.57542E9 FL2 = 1.17645E9 GPRN = GPRN + 400 else: # ERROR!! FL1 = 0 FL2 = 0 GPRN = 0 return FL1, FL2, GPRN def load_gnss(file): with np.load(file, allow_pickle=True) as fid: Gdata = fid['Gdata'] Gdt = fid['Gdt'] Gsats = fid['Gsats'] return Gdata, Gdt, Gsats def read_gnss(dt,prev=False): import os Gdata1 = np.array([]) if prev: for x in ['NRT','RT']: gnssfile = f'../gnss/gnssOrb/GNSS_{dt - timedelta(days=1):%Y.%j}_{x}.npz' if os.path.isfile(gnssfile): if len(Gdata1) == 0: Gdata1, Gdt1, Gsats1 = load_gnss(gnssfile) else: Gdata0, Gdt0, Gsats0 = load_gnss(gnssfile) if Gdt1.shape == Gdt0.shape and np.all(Gdt1 == Gdt0): idx = ~np.isin(Gsats0, Gsats1) Gsats1 = np.hstack([Gsats1, Gsats0[idx]]) Gdata1 = np.hstack([Gdata1, Gdata0[:, idx, :]]) del Gdata0, Gdt0, Gsats0 Gdata2 = np.array([]) for x in ['NRT','RT']: gnssfile = f'../gnss/gnssOrb/GNSS_{dt:%Y.%j}_{x}.npz' if os.path.isfile(gnssfile): if len(Gdata2) == 0: Gdata2, Gdt2, Gsats2 = load_gnss(gnssfile) else: Gdata0, Gdt0, Gsats0 = load_gnss(gnssfile) if Gdt2.shape == Gdt0.shape and np.all(Gdt2 == Gdt0): idx = ~np.isin(Gsats0, Gsats2) Gsats2 = np.hstack([Gsats2, Gsats0[idx]]) Gdata2 = np.hstack([Gdata2, Gdata0[:, idx, :]]) del Gdata0, Gdt0, Gsats0 Gsats = Gsats2 if prev: idx = Gdt1dt_rng[1]) if np.any(tidx): gfiles = gfiles[tidx] tstart = tstart[tidx] tidx = tstart == max(tstart) gfiles = gfiles[tidx] tstart = tstart[tidx] # tend = tend[tidx] Gdt = None Gsats = [] Gdata = [] for ss in src_types: sidx = np.where([ss in gf for gf in gfiles])[0] if sidx.size == 0: continue gfile = gfiles[sidx[0]] Gdata0, Gdt0, Gsats0 = load_gnss(gfile) tidx = np.logical_and(Gdt0>=dt_rng[0]-timedelta(minutes=10),Gdt0 dt_rng[0]) gfiles = gfiles[tidx] tstart = tstart[tidx] Gdt = None Gsats = [] Gdata = [] for gfile in gfiles: Gdata0, Gdt0, Gsats0 = load_gnss(gfile) tidx = np.logical_and(Gdt0>=dt_rng[0]-timedelta(minutes=10),Gdt0dt_rng[1]) gfiles = gfiles[tidx] tstart = tstart[tidx] # tend = tend[tidx] Gdt = None for gsat in gsats: gfile = None for ss in src_types: sidx = np.where([ss in gf for gf in gfiles])[0] if sidx.size==0: continue sfile = gfiles[sidx[np.argmin(np.abs(dt_rng[0] - tstart[sidx] - timedelta(hours=6)))]] with np.load(sfile, allow_pickle=True) as fid: if gsat in fid['Gsats']: gfile = sfile break if gfile is None: raise FileNotFoundError(f'Could not find a GNSS Orb file containing {gsat}.') Gdata0, Gdt0, Gsats0 = load_gnss(gfile) tidx = np.logical_and(Gdt0>=dt_rng[0]-timedelta(minutes=10),Gdt0 proj_lg2[1] exph -= exph[0 if ocd else -1] else: exph = np.zeros_like(l1) range_leo = np.sum(pos_leo_intp ** 2, axis=1) ** 0.5 range_gnss = np.sum(pos_gnss_intp ** 2, axis=1) ** 0.5 delta_rho = GMoCC * np.log((range_gnss + range_leo + range_lg) / (range_gnss + range_leo - range_lg)) tao_old = tao_new.copy() tao_new = (range_lg + exph + delta_rho) / C pos_gnss_intp = lagrange(sec_gnss - pos_gnss[:, 3] / 1E6, pos_gnss[:, :3], sec_ro - tao_new, 9) if not np.any(np.abs(tao_new - tao_old) > 1E-12): break vel_gnss_intp = lagrange(sec_gnss - pos_gnss[:, 3] / 1E6, vel_gnss[:, :3], sec_ro - tao_new, 9) range_dT = np.sum((pos_leo_intp - pos_gnss_intp) ** 2, axis=1) ** 0.5 # from scipy.interpolate import make_interp_spline # speed_dT = make_interp_spline(sec_ro, range_dT, k=2).derivative()(sec_ro) obj_gnss['pos'] = pos_gnss_intp obj_gnss['vel'] = vel_gnss_intp coef_leo = 1 - (GMoCC / range_leo - np.sum(vel_leo_intp ** 2, axis=1) / C / C) / 2 - const propertimedelay_gnss = (coef_leo - 1) * obj_gnss['clk_s'] * C delt_t_rel_gnss = -2 * np.sum(pos_gnss_intp * vel_gnss_intp, axis=1) / C / C delt_t_rel_leo = -2 * np.sum(pos_leo_intp * vel_leo_intp, axis=1) / C / C obj_gnss['range_diff_T'] = range_dT # obj_gnss['range'] = range_lg obj_gnss['GL_term'] = -delta_rho + C * delt_t_rel_gnss + propertimedelay_gnss + obj_gnss['clk_s'] * C obj_leo['delt_t_rel_r'] = delt_t_rel_leo obj_gnss['delt_t_rel_s'] = delt_t_rel_gnss obj_leo['ocd'] = ocd return obj_gnss, obj_leo def remove_cycleslip(ph, ph_m, freq, sat_l, ocd): c0 = 299792458 if not ocd: ph = ph[::-1] ph_m = ph_m[::-1] ph_new = (ph - ph_m) * (2 * freq / c0) if sat_l == 'G1': idx = np.unique(np.where(np.round(np.diff(ph_new[:100])) % 2 == 1)[0] % 2) if idx.size != 1: ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new))))) else: for ii in range(idx[0] + 1, len(ph_new), 2): ph_new[ii:ii + 2] -= np.round(ph_new[ii] - ph_new[ii - 1]) ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new) / 2) * 2))) elif sat_l[0] == 'R': ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new))))) else: ph_new -= np.cumsum(np.concatenate(([0], np.round(np.diff(ph_new) / 2) * 2))) ph_new = ph_new * (c0 / (2 * freq)) + ph_m if not ocd: ph_new = ph_new[::-1] return ph_new # def load_dcb(dt): # import os # from re import sub # import pickle # dcbdir = '../gnss/gnssDCB/' # dcbfile = dcbdir+f'{dt:%Y-%m}_DCB.pkl' # if os.path.exists(dcbfile): # with open(dcbfile, 'rb') as fid: # dcb = pickle.load(fid) # else: # from subprocess import run # from re import search # dcb = {'P1C1': {}, 'P2C2': {}, 'P1P2': {}} # if datetime.now() - dt < timedelta(days=32) and datetime.now().month == dt.month: # zfiles = ["P1C1_RINEX.DCB", # "P2C2_RINEX.DCB", # "P1P2_ALL.DCB"] # rmtdir = "http://ftp.aiub.unibe.ch/CODE/" # else: # zfiles = [f"P1C1{dt:%y%m}_RINEX.DCB.Z", # f"P2C2{dt:%y%m}_RINEX.DCB.Z", # f"P1P2{dt:%y%m}_ALL.DCB.Z"] # rmtdir = f"http://ftp.aiub.unibe.ch/CODE/{dt:%Y}/" # for zfile in zfiles: # run(["wget", rmtdir + zfile], cwd=dcbdir, capture_output=True) # tfile = dcbdir + sub(r'\.Z$', r'', zfile) # if zfile[-1]=='Z': # run(["gunzip", zfile], cwd=dcbdir, capture_output=True) # for line in open(tfile, 'r'): # if len(line) > 47 and search('[GR][0-3][0-9]', line[:3]): # dcb[zfile[:4]][line[:3]] = [float(line[ll:ll + 9].strip()) for ll in [26, 38]] # os.remove(tfile) # with open(dcbfile, 'wb') as fid: # pickle.dump(dcb, fid) # return dcb def calculate_tec(reffile, tecfile=None, dumpnc=True): # import os from re import sub from netCDF4 import Dataset from astropy.coordinates import EarthLocation from zoneinfo import ZoneInfo from subprocess import run # if dumpnc: # out_file = sub(r'(.*)/p01_podPair/(.*)/occTab(.*)\.npz',r'\1/p05_podTec/\2/podTec\3.nc',reffile) # if os.path.exists(out_file) and datetime.fromtimestamp(os.path.getmtime(out_file))>datetime(2026,4,1): # print(f'File {out_file} already exists, skipping.') # return 0,1,1,1 from scipy.interpolate import CubicSpline from frame_conversions import ro_ecef2eci, ro_eci2ecef # from tools import dt2float, get_sat_freq, read_gnss_v2, conv_gnss_ant, phase_correction # Constants c0 = 299792458 kappa = 40.30819294518825 # e^2/(8 pi^2 me eps0) [m^3/s^2] dt0_gps = datetime(1980,1,6) antennas = ['', 'POD1', 'POD2', 'OCC3', 'OCC4'] # Choose reference file dt_ref = [datetime.strptime(x,'%Y%m%d%H%M%S') for x in sub(r'.*(\d{14}\.\d{14}).*',r'\1',reffile.name).split('.')] dt_rr = datetime.strptime(sub(r'.*\.(\d{4}\.\d{3})\..*',r'\1',reffile.name),'%Y.%j') sec_rr = (dt_rr - dt0_gps) / timedelta(seconds=1) # with np.load(reffile, allow_pickle=True) as fid: # data_ref = fid['data'] # yy = data_ref[:, 3] - data_ref[:, 2] # if max(yy)-min(yy)>1000: # raise ValueError('L1 L2 difference out of range!') # Extract info from file names leoName = sub(r'.*occTab_(....)\..*',r'\1',reffile.name) ant_leo_ref = int(sub(r'.*\.A(\d\d)\..*',r'\1',reffile.name)) gnss_ref = sub(r'.*\.(...).npz', r'\1', reffile.name) f1_ref, f2_ref, gprn_ref = get_sat_freq(gnss_ref) c2_ref = 1 / ((f1_ref / f2_ref) ** 2 - 1) c1_ref = c2_ref + 1 filestamp = f'*{leoName}_{dt_ref[0]:%Y%m%d%H%M%S}_{dt_ref[-1]:%Y%m%d%H%M%S}_{gnss_ref}_A{ant_leo_ref:02d}_c*' existing = list(Path('.').glob(f'p05_*/{dt_rr:%Y-%m-%d}/{filestamp}')) if len(existing)>0: raise FileExistsError(f'File {existing[0]} already exists') # ant_gps = np.loadtxt("GPS_Antenna_Infos.Final") # ant_gps_start_time = [datetime(*map(int, x[1:7])) for x in ant_gps] # ant_gps_end_time = [datetime(*map(int, x[7:13])) for x in ant_gps] with np.load('../gnss/GPS_Antenna_Infos.npz', allow_pickle=True) as fid: antoffset_gnss_ref = fid[gnss_ref][()] antoffset_gnss_ref = antoffset_gnss_ref['antenna'][np.where(np.logical_and( antoffset_gnss_ref['dt_rng'][:, 0] < dt_rr, antoffset_gnss_ref['dt_rng'][:, 1] > dt_rr))[0][0]][:3] / 1E3 # Read in gnss position/velocity file accordingly # Gdata, Gdt, Gsats = read_gnss(dt_ref[1] + timedelta(minutes=10), # (dt_ref[0] - timedelta(minutes=10)).day != (dt_ref[1] + timedelta(minutes=10)).day) # if gnss_ro not in Gsats: # raise LookupError(f'RO satellite {gnss_ro} not in Gsats') # if gnss_ref not in Gsats: # raise LookupError(f'REF satellite {gnss_ref} not in Gsats') # gidx = np.hstack([np.where(Gsats == gnss_ro), np.where(Gsats == gnss_ref)])[0] # Gdata = Gdata[:,gidx,:] try: Gdata, Gdt = read_gnss_v2(dt_ref, [gnss_ref]) except: raise LookupError(f'GNSS not found: {gnss_ref}!') tidx = np.logical_and(Gdt > dt_ref[0] - timedelta(minutes=10), Gdt < dt_ref[1] + timedelta(minutes=10)) dts_gnss = Gdt[tidx] sec_gnss = dt2float(dts_gnss, dt_rr) pos_gnss = Gdata[tidx,:,:4] vel_gnss = Gdata[tidx,:,4:] del Gdata, Gdt # Read in leo position/velocity/attitude files accordingly with np.load(f"p01_leoOrb/leoOrb_{leoName}.{dt_rr:%Y.%j}.npz", allow_pickle=True) as fid: Ldata = fid['data'] Ldts = fid['dts'] tidx = np.logical_and(Ldts > dt_ref[0] - timedelta(minutes=2), Ldts < dt_ref[1] + timedelta(minutes=2)) dts_leoorb = Ldts[tidx] pos_leo = Ldata[tidx, :4] vel_leo = Ldata[tidx, 4:] # tmsk = ~np.isin(np.arange(len(Ldts)), # np.unique(np.where(np.diff(Ldts) > timedelta(seconds=5))[0] + np.arange(-30, 90)[:, None])) # if np.all(Ldts[tmsk] < dt_ref[1]) or np.all(Ldts[tmsk] > dt_ref[0]): # raise ValueError( # f'LEO time ({Ldts[tmsk][0]:%d-%H:%M} >> {Ldts[tmsk][-1]:%d-%H:%M}) could not cover REF event ({dt_ref[0]:%H:%M} >> {dt_ref[1]:%H:%M})!') # tidx = np.logical_and(Ldts > dt_ref[0] - timedelta(minutes=1), Ldts < dt_ref[1] + timedelta(minutes=1)) # if np.sum(tidx)<20: # raise ValueError( # f'LEO with large gap, no enough point in REF event ({dt_ref[0]:%H:%M} >> {dt_ref[1]:%H:%M})!') # if np.all(tmsk[tidx]): # dts_leoorb = Ldts[tidx] # pos_leo = Ldata[tidx,:4] # vel_leo = Ldata[tidx,4:] # else: # tidx0 = np.where(np.logical_and(Ldts < dt_ref[0], tmsk))[0] # tidx1 = np.where(np.logical_and(Ldts > dt_ref[1], tmsk))[0] # tidx = [tidx0[-60 if len(tidx0)>60 else 0],tidx1[min(60, len(tidx1)-1)]] # dts_leoorb = Ldts[tidx[0]:tidx[1]] # pos_leo = Ldata[tidx[0]:tidx[1], :4] # vel_leo = Ldata[tidx[0]:tidx[1], 4:] # tmsk = tmsk[tidx[0]:tidx[1]] # sec_leoorb = dt2float(dts_leoorb, dt_rr) # sec_mid = (sec_leoorb[0] + sec_leoorb[-1]) / 2 # sec_rng = (sec_leoorb[-1] - sec_leoorb[0]) / 2 # dts_leoorb = np.arange(dts_leoorb[0],dts_leoorb[-1],timedelta(seconds=1)).astype(datetime) # pos_leo = np.array([np.polyval(np.polyfit((sec_leoorb[tmsk] - sec_mid) / sec_rng, x, 10), # (np.arange(sec_leoorb[0], sec_leoorb[-1]) - sec_mid) / sec_rng) for x in pos_leo[tmsk, :3].T]+ # [np.interp(np.arange(sec_leoorb[0], sec_leoorb[-1]), sec_leoorb, pos_leo[:, 3])]).T # vel_leo = np.array([np.polyval(np.polyfit((sec_leoorb[tmsk] - sec_mid) / sec_rng, x, 10), # (np.arange(sec_leoorb[0], sec_leoorb[-1]) - sec_mid) / sec_rng) for x in vel_leo[tmsk, :3].T]+ # [np.interp(np.arange(sec_leoorb[0], sec_leoorb[-1]), sec_leoorb, vel_leo[:, 3])]).T del Ldata, Ldts with np.load(f"p01_leoAtt/leoAtt_{leoName}.{dt_rr:%Y.%j}.npz", allow_pickle=True) as fid: Ldata = fid['data'] Ldts = fid['dts'] Lvar = fid['vars'] Ldts, tidx = np.unique(Ldts, return_index=True) Ldata = Ldata[tidx] tidx = np.logical_and(Ldts > dt_ref[0] - timedelta(minutes=1), Ldts < dt_ref[1] + timedelta(minutes=1)) if ~np.any(tidx): raise ValueError('No matching leoOrb/leoAtt') dts_leoatt = Ldts[tidx] att_leo = Ldata[tidx, :][:, ['att' in x for x in Lvar]] sad_leo = Ldata[tidx, :][:, ['sad' in x for x in Lvar]] aridx = np.cumsum(np.hstack(([False], np.sum(att_leo[1:] * att_leo[:-1], axis=1) < 0))) % 2 == 1 att_leo[aridx] *= -1 del Ldata, Ldts, Lvar tidx = np.logical_and(dts_leoorb >= dts_leoatt[0], dts_leoorb <= dts_leoatt[-1]) dts_leo = dts_leoorb[tidx] sec_leo = dt2float(dts_leo, dt_rr) pos_leo = pos_leo[tidx] vel_leo = vel_leo[tidx] _, idx = np.unique(dts_leoatt, return_index=True) att_leo = CubicSpline(dt2float(dts_leoatt[idx],dts_leo[0]), att_leo[idx])(dt2float(dts_leo,dts_leo[0])) att_leo /= np.sum(att_leo ** 2, axis=1)[:,None] ** 0.5 sad_leo -= np.cumsum(np.round(np.diff(sad_leo, axis=0, prepend=sad_leo[:1]) / 360) * 360, axis=0) sad_leo = CubicSpline(dt2float(dts_leoatt[idx],dts_leo[0]), sad_leo[idx])(dt2float(dts_leo,dts_leo[0])) sad_leo[:, 0] = sad_leo[:, 0] % 360 sad_leo[:, 1] = (sad_leo[:, 1] + 180) % 360 - 180 del dts_leoorb, dts_leoatt # tsmp = 1 # 30 # tspc = np.round(np.linspace(0,len(sec_leo)-1,min(max(-((len(sec_leo)-1)//-tsmp)+1,20),len(sec_leo)))).astype(int) # dts_leo = dts_leo[tspc] # sec_leo = sec_leo[tspc] # pos_leo = pos_leo[tspc] # vel_leo = vel_leo[tspc] # att_leo = att_leo[tspc] # Load ro/ref files with np.load(reffile, allow_pickle=True) as fid: data_ref = fid['data'] header_ref = np.hstack([fid['header'],['C1C','C2X']]) l_ref = np.where([x[0] == 'L' for x in header_ref])[0] s_ref = np.where([x[0] == 'S' for x in header_ref])[0] c_ref = np.where([x[0] == 'C' for x in header_ref])[0] if np.any(np.abs(data_ref[:, c_ref[1]] - data_ref[:, c_ref[0]]) > 100): raise ValueError( f"Out-of-range pseudorange difference detected: Min={np.min(data_ref[:, c_ref[1]] - data_ref[:, c_ref[0]]):.2f} Max={np.max(data_ref[:, c_ref[1]] - data_ref[:, c_ref[0]]):.2f}") idx_seg = np.where(np.abs(np.diff(data_ref[:, l_ref[1]] - data_ref[:, l_ref[0]])) > 0.5)[0] if len(idx_seg) > 0: idx_seg = np.r_[0, idx_seg + 1, data_ref.shape[0]] idx_long = np.argmax(np.diff(idx_seg)) data_ref = data_ref[idx_seg[idx_long]:idx_seg[idx_long + 1]] dt_ref = [dt0_gps + timedelta(seconds=data_ref[0, 0] + data_ref[0, 1]), dt0_gps + timedelta(seconds=data_ref[-1, 0] + data_ref[-1, 1])] reffile.unlink() reffile = reffile.with_stem( f"occTab_{leoName}.{dt_rr:%Y.%j}.{dt_ref[0]:%Y%m%d%H%M%S}.{dt_ref[1]:%Y%m%d%H%M%S}.A{ant_leo_ref:02d}.{gnss_ref}") np.savez(reffile, data=data_ref, header=header_ref) sec_ref = data_ref[:, 0] - sec_rr + data_ref[:, 1] if sec_leo[0]>sec_ref[0] or sec_leo[-1]> {sec_leo[-1]:.1f}) could not cover REF event ({sec_ref[0]:.1f} >> {sec_ref[-1]:.1f})!') sec_ref -= CubicSpline(sec_leo, pos_leo[:,-1])(sec_ref)/1E6 if ~np.all(np.isfinite(pos_gnss[:, 0, 3])): raise ValueError(f'{np.sum(~np.isfinite(pos_gnss[:, 0, 3]))}/{sec_gnss.size} clock bias value for {gnss_ref} is not finite.') ant_ref_ecef = conv_gnss_ant(pos_gnss[:, 0, :3], antoffset_gnss_ref, dts_gnss) pos_leo_eci, vel_leo_eci = ro_ecef2eci(pos_leo, vel_leo, dts_leo) pos_ref_eci, vel_ref_eci = ro_ecef2eci(pos_gnss[:, 0, :3] + ant_ref_ecef, vel_gnss[:, 0, :], dts_gnss) pos_leo_eci = np.hstack([pos_leo_eci, pos_leo[:,[-1]]]) vel_leo_eci = np.hstack([vel_leo_eci, pos_leo[:, [-1]]]) pos_ref_eci = np.hstack([pos_ref_eci, pos_gnss[:, 0, [3]]]) vel_ref_eci = np.hstack([vel_ref_eci, vel_gnss[:, 0, [3]]]) with np.load('antenna_offset.npz') as fid: antoffset_leo_ref = fid[f'ANT_{leoName}_{antennas[ant_leo_ref]}'][12:15] antbsight_leo = fid[f'ANT_{leoName}_{antennas[ant_leo_ref]}'][15:18] obj_gnss_ref, obj_leo_ref = phase_correction(data_ref[:, l_ref[0]], pos_ref_eci, vel_ref_eci, sec_gnss, pos_leo_eci, vel_leo_eci, att_leo, sec_leo, sec_ref, antoffset_leo_ref, True) phL1_ref = data_ref[:, l_ref[0]] - obj_gnss_ref['range_diff_T'] + obj_gnss_ref['GL_term'] # coef = np.polyfit(sec_ref, phL1_ref, 1) # if leoName.startswith('YM') and (np.abs(coef[0]) > 1 or np.max(np.abs(phL1_ref - np.polyval(coef, sec_ref))) > 10): # raise ValueError('L1 non-linear!') phL2_ref = data_ref[:, l_ref[1]] - obj_gnss_ref['range_diff_T'] + obj_gnss_ref['GL_term'] dph_ref = phL1_ref - phL2_ref # dph_ref = np.hstack( # [np.cumsum(dph_ref[:49])[::2] / np.arange(1, 51, 2), np.convolve(dph_ref, np.ones((51,)), 'valid') / 51, # (np.cumsum(dph_ref[:-50:-1])[::2] / np.arange(1, 51, 2))[::-1]]) termL1_ref = phL1_ref + c2_ref * dph_ref termL2_ref = phL2_ref + c1_ref * dph_ref # termL1_ro = obj_leo_ro['clk_m'] + obj_leo_ro['delt_t_rel_r'] * c0 # Alternative method # termL2_ro = obj_leo_ro['clk_m'] + obj_leo_ro['delt_t_rel_r'] * c0 exL1 = phL1_ref # - termL1_ref exL2 = phL2_ref # - termL2_ref exL1 -= exL1[0 if obj_leo_ref['ocd'] else -1] exL2 -= exL2[0 if obj_leo_ref['ocd'] else -1] OL = obj_leo_ref['pos'] OG = obj_gnss_ref['pos'] if np.all(np.isnan(OG)) or np.all(np.isnan(OL)): raise ValueError('No finite coordinate values') LG = OG - OL antbsight_eci = inv_quat_rotate(att_leo, antbsight_leo) ex = lagrange(sec_leo / 86400, antbsight_eci, np.round(sec_ref) / 86400, 11) ex /= np.sum(ex ** 2, axis=1)[:, None] ** 0.5 ez = -OL # * [1, 1, 1.00669437999014] ey = np.cross(ez, ex) ey /= np.sum(ey ** 2, axis=1)[:, None] ** 0.5 ez = np.cross(ex, ey) LG_bsight = np.array([np.sum(LG * x, axis=1) for x in [ex, ey, ez]]).T el_bsight = -np.arcsin(LG_bsight[:, 2] / np.sum(LG_bsight ** 2, axis=1) ** 0.5) / np.pi * 180 az_bsight = np.arctan2(LG_bsight[:, 1], LG_bsight[:, 0]) / np.pi * 180 # sad_ref = lagrange(sec_leo / 86400, sad_leo, np.round(sec_ref) / 86400, 11) sad_leo -= np.cumsum(np.round(np.diff(sad_leo, axis=0, prepend=sad_leo[:1]) / 360) * 360, axis=0) sad_ref = lagrange(dt2float(dts_leo, dt_rr) / 86400, sad_leo, sec_ref / 86400, 1) sad_ref[:, 0] = sad_ref[:, 0] % 360 sad_ref[:, 1] = (sad_ref[:, 1] + 180) % 360 - 180 if gnss_ref[0] == 'G': bslut = list((working_dir / 'p03_bslut').glob(f'bslut_{leoName}.*.D*.A{ant_leo_ref:02d}.GL1.npz'))[0] with np.load(bslut, allow_pickle=True) as fid: mpe_grid = fid['mpe_grid'][:,:,:,:,:1] import pickle with open('p01_podCode/podCode_Psats.pkl','rb') as fid: ptimes = pickle.load(fid) leoidx = int(leoName[3])-1 if gnss_ref in ptimes and (ptimes[gnss_ref][leoidx] == None or dt_ref[0] < ptimes[gnss_ref][leoidx]): bslut = list((working_dir / 'p03_bslut').glob(f'bslut_{leoName}.*.D*.A{ant_leo_ref:02d}.GP2.npz'))[0] else: bslut = list((working_dir / 'p03_bslut').glob(f'bslut_{leoName}.*.D*.A{ant_leo_ref:02d}.GC2.npz'))[0] with np.load(bslut, allow_pickle=True) as fid: mpe_grid = np.concatenate((mpe_grid,fid['mpe_grid'][:,:,:,:,1:]),axis=4) elif gnss_ref[0] == 'R': gnssidx = int(gnss_ref[1:]) gnssidx -= 4 if (gnssidx - 1) % 8 > 3 else 0 bslut = list((working_dir / 'p03_bslut').glob(f'bslut_{leoName}.*.D*.A{ant_leo_ref:02d}.R{gnssidx:02d}.npz'))[0] with np.load(bslut, allow_pickle=True) as fid: mpe_grid = fid['mpe_grid'] if gnss_ref[0] in 'GR': sda_bins = np.arange(-180,180,10) mpe_grid = mpe_grid[((el_bsight + 90)//2).astype('int'), ((az_bsight + 180)//2%180).astype('int'), ((sad_ref[:,1]+45)//30).astype('int')] idx = np.any(np.isfinite(mpe_grid), axis=(0, 2)) mpe = np.array( [[np.interp(sad_ref[x,0], sda_bins[idx], mpe_grid[x, idx, y], period=360) for y in range(2)] for x in range(len(el_bsight))]) # mpe[mpe>0.3] = 0.3 # mpe[mpe<-0.3] = -0.3 else: mpe = np.zeros((len(el_bsight),2)) # TC: updated 2026/06/04 (4 lines, get the coefficient for vector addition, remove the extrapolation in altp_km) pcoef = np.sum(OL * LG, axis=1) / np.sum(LG * LG, axis=1) proj_lg = OL - LG * pcoef[:, None] altp_km = EarthLocation.from_geocentric(*proj_lg.T, unit='m').geodetic.height.value / 1E3 altp_km[pcoef > 0] = np.nan altl_km = EarthLocation.from_geocentric(*OL.T, unit='m').geodetic.height.value / 1E3 newZ = proj_lg[0,:].copy() newY = np.cross(OL[0,:],OG[0,:]) newX = np.cross(newY,newZ) newX /= sum(newX ** 2) ** .5 newY /= sum(newY ** 2) ** .5 newZ /= sum(newZ ** 2) ** .5 mrot = np.array([newX,newY,newZ]) # cosPL = np.sum(proj_lg*OL,axis=1)/np.sum(OL*OL,axis=1)**.5/np.sum(proj_lg*proj_lg,axis=1)**.5 # cosPG = np.sum(proj_lg*OG,axis=1)/np.sum(OG*OG,axis=1)**.5/np.sum(proj_lg*proj_lg,axis=1)**.5 # RL = (np.sum(OL ** 2, axis=1) ** .5)[:, None] * np.array([-(1 - cosPL ** 2) ** .5, cosPL]).T / 6378131 # RG = (np.sum(OG ** 2, axis=1) ** .5)[:, None] * np.array([ (1 - cosPG ** 2) ** .5, cosPG]).T / 6378131 RL = (mrot@OL.T).T/6378137 RG = (mrot@OG.T).T/6378137 coszn = np.sum(OL * LG, axis=1) / np.sum(LG * LG, axis=1)**0.5 / np.sum(OL * OL, axis=1)**0.5 # if np.all(coszn<0): # raise ValueError('All observation below horizon!') alti_km = 4000 # km rad_ratio = (alti_km+6371)/(altl_km+6371) meps = (coszn + (rad_ratio**2+coszn**2-1)**0.5)/(1+rad_ratio) # idx = np.where(np.isfinite(exL1))[0] # exL1 -= exL1[idx[0 if obj_leo_ref['ocd'] else -1]] # exL2 -= exL2[idx[0 if obj_leo_ref['ocd'] else -1]] # exLC = exL1 + c2_ref * (exL1 - exL2) tt = np.arange(mpe.shape[0]) idx = np.all(np.isfinite(mpe), axis=1) if np.any(idx): mpe = np.array([np.interp(tt,tt[idx],mpe[idx,x]) for x in [0,1]]).T else: mpe[:, :] = 0 Pdiff = data_ref[:,l_ref[0]] - data_ref[:,l_ref[1]] fr_ref = 2 / ((f1_ref / f2_ref) ** 2 - 1) + 1 S = data_ref[:, s_ref] ** 2 # S[coszn<0] = 0 D1 = data_ref[:, c_ref[0]] - fr_ref * data_ref[:, l_ref[0]] + (fr_ref - 1) * data_ref[:, l_ref[1]] - mpe[:,0] D2 = data_ref[:, c_ref[1]] - (fr_ref + 1) * data_ref[:, l_ref[0]] + fr_ref * data_ref[:, l_ref[1]] - mpe[:,1] S /= np.sum(S, axis=0)[None, :] Ddiff = np.nansum(D1 * S[:, 0]) - np.nansum(D2 * S[:, 1]) Cdiff = (data_ref[:, c_ref[1]] - data_ref[:, c_ref[0]]) / (f2_ref ** -2 - f1_ref ** -2) / kappa / 1E16 # tec = (Pdiff - Ddiff) / (f2_ref ** -2 - f1_ref ** -2) / kappa / 1E16 tec = (Pdiff - Ddiff) / (f2_ref ** -2 - f1_ref ** -2) / kappa / 1E16 # tec0 = tec.copy() tec -= np.nanmean(tec-Cdiff) leveling_offset = -Ddiff / (f2_ref ** -2 - f1_ref ** -2) / kappa / 1E16 leveling_err = np.sum((tec - Cdiff) ** 2) ** 0.5 / len(tec) # dcb = load_dcb(dt_ref[0]) # dcb_factor = c0 / (f2_ref ** -2 - f1_ref ** -2) / kappa / 1E16 # dcb_value = (dcb['P1P2'][gnss_ref][0] - dcb['P1C1'][gnss_ref][0] + dcb['P2C2'][gnss_ref][0]) / 1E9 * dcb_factor # tec += dcb_value dts_ref = [dt_rr + timedelta(seconds=x) for x in sec_ref] pos_leo_ecef, vel_leo_ecef, pos_leo_tod, vel_leo_tod = ro_eci2ecef(obj_leo_ref['pos'], obj_leo_ref['vel'], dts_ref, tod=True) pos_gnss_ecef, vel_gnss_ecef, pos_gnss_tod, vel_gnss_tod = ro_eci2ecef(obj_gnss_ref['pos'], obj_gnss_ref['vel'], dts_ref, tod=True) with np.load(reffile, allow_pickle=True) as fid: data = fid['data'] header = fid['header'] np.savez(reffile, data=data, header=header, sec_ref=sec_ref, tec=tec, coszn=coszn, meps=meps) if dumpnc: if leveling_err > 1: raise ValueError(f'Leveling error is too large: {leveling_err}') # TC: updated 2026/06/01 (removed 3 lines) # with np.load('dcb_default.npz', allow_pickle=True) as fid: # allgnsses = fid['allgnsses'][0 if leoName == 'GN04' else 1, ant_leo_ref] # alldcbs = fid['alldcbs'][0 if leoName == 'GN04' else 1, ant_leo_ref] # with np.load('dcb_daily.npz', allow_pickle=True) as fid: # allgnsses = fid['allgnsses'] # alldcbs = fid['alldcbs'][0 if leoName == 'GN04' else 1, ant_leo_ref] # alldts = fid['alldts'] # alldcbs = alldcbs[dt_rr==alldts][0] #TC: updated 2026/06/01 (1 line, change default value) dcb = np.nan for dayprev in range(1,366): dt_now = dt_rr - timedelta(days=dayprev) dcbfile = Path(f'p04_teczdcb/teczdcb_{leoName}.{dt_now:%Y-%m-%d}.A{ant_leo_ref:02d}.npz') if dcbfile.is_file(): with np.load(dcbfile, allow_pickle=True) as fid: allgnsses = fid['gnsses'] alldcbs = fid['dcb_sol'] if gnss_ref in allgnsses: dcb = alldcbs[allgnsses==gnss_ref][0] break if np.any(tec<-100-dcb): raise ValueError(f'TEC too negative: Min={np.nanmin(tec)+dcb} Max={np.nanmax(tec)+dcb}') if np.any(tec>1000-dcb): raise ValueError(f'TEC too large: Min={np.nanmin(tec)+dcb} Max={np.nanmax(tec)+dcb}') # TC: updated 2026-06-01 (added 8 lines, check existence of dcb value) if np.isfinite(dcb): out_dir = working_dir / f'p05_podTec/{dt_rr:%Y-%m-%d}' edp_dir = working_dir / f'p06_podEdp/{dt_rr:%Y-%m-%d}' edp_dir.mkdir(exist_ok=True,parents=True) else: out_dir = working_dir / f'p05_forDcb/{dt_rr:%Y-%m-%d}' dcb = 0.0 out_dir.mkdir(exist_ok=True) filestamp = f'{leoName}_{dt_ref[0]:%Y%m%d%H%M%S}_{dt_ref[-1]:%Y%m%d%H%M%S}_{gnss_ref}_A{ant_leo_ref:02d}_c{datetime.now().astimezone(ZoneInfo("UTC")):%Y%m%d%H%M%S}' out_file = out_dir / f'ROSW_TEC_{filestamp}_nes_ops.nc' dt_now = datetime.now() # TC: updated 2026/06/04 (remove 5 lines) idx = [np.argmax(coszn), np.argmax(tec), 0, -1] pos_attr = EarthLocation.from_geocentric(*pos_leo_ecef[idx].T, unit='m').geodetic alt_attr = pos_attr.height.value / 1E3 lat_attr = pos_attr.lat.value lon_attr = pos_attr.lon.value lct_attr = (lon_attr/15 + sec_ref[idx]/3600) % 24 with Dataset(out_file, 'w') as fid: fid.createDimension('time',len(tec)) fid.createDimension('attribute_dimension',1) fid.createVariable('time', 'f8', ('time', )) fid['time'].setncatts( {'_FillValue': -999.0, 'calendar': 'standard', 'long_name': f"Seconds since start time {sec_rr + sec_ref[0]:.1f}", 'units': f"Seconds since start time {sec_rr + sec_ref[0]:.1f}", 'add_offset': sec_rr + sec_ref[0]} ) fid['time'][...] = sec_ref+sec_rr fid.createVariable('TEC', 'f8', ('time',)) fid['TEC'].setncatts( {'_FillValue': -999.0, 'units': 'TECU', 'long_name': 'Total Electron Content along LEO-GPS link'} ) fid['TEC'][...] = tec+dcb fid.createVariable('antid', 'int16', ('time',)) fid['antid'].setncatts( {'_FillValue': np.int16(-999), 'units': 'Unitless', 'long_name': 'Antenna ID for each time step'} ) fid['antid'][...] = np.full((len(tec),),ant_leo_ref).astype('int16') fid.createVariable('elevation', 'f8', ('time',)) fid['elevation'].setncatts( {'_FillValue': -999.0, 'units': 'deg', 'long_name': 'Elevation angle of LEO-GPS link'} ) fid['elevation'][...] = np.arcsin(coszn)/np.pi*180 fid.createVariable('occheight', 'f8', ('time',)) fid['occheight'].setncatts( {'_FillValue': -999.0, 'units': 'km', 'long_name': 'Occultation perigee point height'} ) fid['occheight'][...] = altp_km fid.createVariable('caL1_SNR', 'f4', ('time',)) fid['caL1_SNR'].setncatts( {'_FillValue': np.float32(-999.0), 'long_name': 'First frequency signal strength'} ) fid['caL1_SNR'][...] = data_ref[:, s_ref[0]].astype('float32') fid.createVariable('pL2_SNR', 'f4', ('time',)) fid['pL2_SNR'].setncatts( {'_FillValue': np.float32(-999.0), 'long_name': 'Second frequency signal strength'} ) fid['pL2_SNR'][...] = data_ref[:, s_ref[1]].astype('float32') fid.createVariable('x_LEO', 'f8', ('time',)) fid['x_LEO'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF LEO satellite X position', 'valid_range': [-8000.0, 8000.0]} ) fid['x_LEO'][...] = pos_leo_ecef[:, 0] / 1E3 fid.createVariable('y_LEO', 'f8', ('time',)) fid['y_LEO'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF LEO satellite Y position', 'valid_range': [-8000.0, 8000.0]} ) fid['y_LEO'][...] = pos_leo_ecef[:, 1] / 1E3 fid.createVariable('z_LEO', 'f8', ('time',)) fid['z_LEO'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF LEO satellite Z position', 'valid_range': [-8000.0, 8000.0]} ) fid['z_LEO'][...] = pos_leo_ecef[:, 2] / 1E3 fid.createVariable('x_GNSS', 'f8', ('time',)) fid['x_GNSS'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF GNSS satellite X position', 'valid_range': [-30000.0, 30000.0]} ) fid['x_GNSS'][...] = pos_gnss_ecef[:, 0] / 1E3 fid.createVariable('y_GNSS', 'f8', ('time',)) fid['y_GNSS'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF GNSS satellite Y position', 'valid_range': [-30000.0, 30000.0]} ) fid['y_GNSS'][...] = pos_gnss_ecef[:, 1] / 1E3 fid.createVariable('z_GNSS', 'f8', ('time',)) fid['z_GNSS'].setncatts( {'_FillValue': -99999.0, 'units': 'km', 'long_name': 'ECEF GNSS satellite Z position', 'valid_range': [-30000.0, 30000.0]} ) fid['z_GNSS'][...] = pos_gnss_ecef[:, 2] / 1E3 fid.createVariable('exph_L1', 'f8', ('time',)) fid['exph_L1'].setncatts( {'units': 'm', 'long_name': 'Excess phase L1'} ) fid['exph_L1'][...] = exL1 fid.createVariable('exph_L2', 'f8', ('time',)) fid['exph_L2'].setncatts( {'units': 'm', 'long_name': 'Excess phase L2'} ) fid['exph_L2'][...] = exL2 fid.createVariable('bsangle', 'f8', ('time',)) fid['bsangle'].setncatts( {'units': 'Degrees', 'long_name': 'Angle between boresight vector and LEO->GNSS vector'} ) fid['bsangle'][...] = np.arccos(LG_bsight[:, 0] / np.sum(LG_bsight ** 2, axis=1) ** 0.5) / np.pi * 180 fid.createVariable('leveling_offset', 'f8', ('attribute_dimension',)) fid['leveling_offset'].setncatts( {'units': 'TEC Units', 'long_name': 'Phase/code leveling offset: Absolute TEC = relative TEC + leveling_offset + GNSS DCB + LEO DCB'} ) fid['leveling_offset'][...] = leveling_offset fid.setncatts({ 'start_time': np.uint32(np.round(sec_ref[0]+sec_rr)), 'stop_time': np.uint32(np.round(sec_ref[-1]+sec_rr)), 'gobs_start': np.uint32(np.round(sec_ref[0]+sec_rr)), 'gobs_stop': np.uint32(np.round(sec_ref[-1]+sec_rr)), 'History': f'Created by atec.py {dt_now:%d/%m/%Y %H:%M:%S}', 'Title': 'Absolute TEC measurements from radio occultation', 'Conventions': 'CF-1.6', 'year': np.int32(dts_ref[0].year), 'month': np.int32(dts_ref[0].month), 'day': np.int32(dts_ref[0].day), 'hour': np.int32(dts_ref[0].hour), 'minute': np.int32(dts_ref[0].minute), 'second': dts_ref[0].second+dts_ref[0].microsecond/1E6, 'tecmin': np.min(tec+dcb), 'tecmax': np.max(tec+dcb), 'elevmin': np.min(np.arcsin(coszn)/np.pi*180), 'elevmax': np.max(np.arcsin(coszn)/np.pi*180), 'occaltmin': np.nanmin(altp_km) if np.any(np.isfinite(altp_km)) else np.int32(-999), 'occaltmax': np.nanmax(altp_km) if np.any(np.isfinite(altp_km)) else np.int32(-999), 'arcTab_start': np.uint32(np.round(sec_ref[0]+sec_rr)), 'arcTab_stop': np.uint32(np.round(sec_ref[-1]+sec_rr)), 'arcTab_ID': 1, 'alt_elevmax': alt_attr[0], 'lat_elevmax': lat_attr[0], 'lon_elevmax': lon_attr[0], 'lct_elevmax': lct_attr[0], 'alt_tecmax': alt_attr[1], 'lat_tecmax': lat_attr[1], 'lon_tecmax': lon_attr[1], 'lct_tecmax': lct_attr[1], 'alt_start': alt_attr[2], 'lat_start': lat_attr[2], 'lon_start': lon_attr[2], 'lct_start': lct_attr[2], 'alt_stop': alt_attr[3], 'lat_stop': lat_attr[3], 'lon_stop': lon_attr[3], 'lct_stop': lct_attr[3], 'processing_center': 'NESDIS', 'L2C': 1, 'creation_time': f'{dt_now:%Y-%m-%d %H:%M:%S.%f}', 'mission': 'gnomestec', 'dump_id': f'{dt_ref[0]:%Y.%j}.0{leoName[2:]}.{ant_leo_ref:02d}', 'leo_id': np.int32(leoName[3]), 'prn_id': np.int32(gnss_ref[1:]), 'duration': sec_ref[-1]-sec_ref[0], 'leodcb_flag': np.int32(1), 'leodcb': dcb, 'gpsdcb_flag': np.int32(0), 'gpsdcb': np.int32(-999), 'leveling_offset': leveling_offset, 'leveling_err': leveling_err, 'conid': gnss_ref[0], 'dcb_units': 'TECU', }) if np.isfinite(dcb): results = run([(working_dir / '.venv/bin/python').as_posix(), 'TEC_to_EDP_final_gnss.py', '--TEC', out_file.as_posix(), '--out', edp_dir.as_posix()], capture_output=True, cwd=working_dir.parent / 'TEC_EDP') return sec_ref, tec, coszn, meps def combine_dcb(): from glob import glob from re import sub gnsses = np.array( ['C06', 'C07', 'C08', 'C09', 'C10', 'C11', 'C12', 'C13', 'C14', 'C16', 'C19', 'C20', 'C21', 'C22', 'C23', 'C24', 'C25', 'C26', 'C27', 'C28', 'C29', 'C30', 'C32', 'C33', 'C34', 'C35', 'C36', 'C37', 'C38', 'C39', 'C40', 'C41', 'C42', 'C43', 'C44', 'C45', 'C47', 'C48', 'C49', 'C50', 'E02', 'E03', 'E04', 'E05', 'E06', 'E07', 'E08', 'E09', 'E10', 'E11', 'E12', 'E13', 'E14', 'E15', 'E16', 'E18', 'E19', 'E21', 'E23', 'E24', 'E25', 'E26', 'E27', 'E29', 'E30', 'E31', 'E33', 'E34', 'E36', 'G01', 'G02', 'G03', 'G04', 'G05', 'G06', 'G07', 'G08', 'G09', 'G10', 'G11', 'G12', 'G13', 'G14', 'G15', 'G16', 'G17', 'G18', 'G19', 'G20', 'G21', 'G22', 'G23', 'G24', 'G25', 'G26', 'G27', 'G28', 'G29', 'G30', 'G31', 'G32', 'J02', 'J03', 'J04', 'R02', 'R03', 'R04', 'R05', 'R06', 'R07', 'R08', 'R09', 'R11', 'R12', 'R14', 'R15', 'R16', 'R17', 'R18', 'R19', 'R20', 'R21', 'R22', 'R24', 'R25', 'R26', 'R27', 'R28']) ants = [0,1] leos = ['GN04','GN05'] # gnsstypes = ['C','E'] alldcbs = [[[],[]],[[],[]]] allgnsses = [[[],[]],[[],[]]] alldts = [[[],[]],[[],[]]] kernal = np.array([1]) # np.arange(10, 0, -1) / 55 for aa,ant in enumerate(ants): for ll,leo in enumerate(leos): dcbfiles = glob(f'p04_teczdcb/teczdcb_{leo}.*.A{ant:02d}.npz') dcbfiles.sort() dts = np.array([datetime.strptime(sub(r'.*\.(\d{4}-\d\d-\d\d)\..*',r'\1',x),'%Y-%m-%d') for x in dcbfiles]) idx = dts>datetime(2025,12,15) dcbfiles = np.array(dcbfiles)[idx] dts = dts[idx] alldcb = np.full((len(dts),len(gnsses)),np.nan) for dd, dt, dcbfile in zip(range(len(dts)),dts,dcbfiles): with np.load(dcbfile, allow_pickle=True) as fid: # sec_tec = fid['sec_tec'] # tecz = fid['tecz'] dcb_gnss = fid['gnsses'] dcb_sol = fid['dcb_sol'] alldcb[dd,np.isin(gnsses,dcb_gnss)] = dcb_sol idx = np.any(np.isfinite(alldcb), axis=0) allgnss = gnsses[idx] alldcb = alldcb[:, idx] idx = np.any(np.abs(alldcb-np.nanmedian(alldcb,axis=0))>3.5,axis=1) alldcb[idx] = np.nan for ii in range(alldcb.shape[1]): idx = np.isfinite(alldcb[:,ii]) alldcb[:,ii] = np.interp(dt2float(dts,dts[0]),dt2float(dts[idx],dts[0]),alldcb[idx,ii]) alldcbs[ll][aa] = np.array([np.convolve(kernal,alldcb[:,x],'valid') for x in range(alldcb.shape[1])]).T # alldcbs[ll][aa] = np.nanmedian(alldcb, axis=0) allgnsses[ll][aa] = allgnss.copy() alldts[ll][aa] = dts[len(kernal)-1:] cmb_gnss = [] cmb_dts = [] for aa,ant in enumerate(ants): for ll,leo in enumerate(leos): cmb_gnss = np.union1d(cmb_gnss,allgnsses[aa][ll]) cmb_dts = np.union1d(cmb_dts,alldts[aa][ll]) cmb_dcbs = np.full((len(ants),len(leos),len(cmb_dts),len(cmb_gnss)),np.nan) for aa,ant in enumerate(ants): for ll,leo in enumerate(leos): tidx = np.isin(cmb_dts,alldts[aa][ll]) gidx = np.isin(cmb_gnss,allgnsses[aa][ll]) idx = np.ix_([aa],[ll],tidx,gidx) cmb_dcbs[idx] = alldcbs[aa][ll] idx = np.isfinite(cmb_dcbs) if ~np.all(idx): for tt,dts in enumerate(cmb_dts): for gg,gnss in enumerate(cmb_gnss): if np.any(~np.isfinite(cmb_dcbs[:,:,tt,gg])): for aa,ant in enumerate(ants): idx = np.where(np.isfinite(cmb_dcbs[aa, :, tt, gg]))[0] if len(idx)>0: cmb_dcbs[aa, :, tt, gg] = np.interp(np.arange(len(leos)), idx, cmb_dcbs[aa, idx, tt, gg]) for ll,leo in enumerate(leos): idx = np.where(np.isfinite(cmb_dcbs[:, ll, tt, gg]))[0] if len(idx)>0: cmb_dcbs[:, ll, tt, gg] = np.interp(np.arange(len(ants)), idx, cmb_dcbs[idx, ll, tt, gg]) np.savez('dcb_daily.npz',alldts=cmb_dts,alldcbs=cmb_dcbs,allgnsses=cmb_gnss) return def apply_dcb(args): iprof = args[0] dt = args[1] leo = args[2] ant = args[3] out_tag = '' if len(args)>4: out_tag = args[4] from netCDF4 import Dataset from re import sub from shutil import copy # dts = [datetime(2026, 1, 1) + timedelta(days=dd) for dd in range(4)] # leos = ['GN04', 'GN05'] # ants = [0, 1] # for dt in dts: # for leo in leos: # for ant in ants: dcbfile = working_dir / f'p04_teczdcb/teczdcb_{leo}.{dt-timedelta(days=1):%Y-%m-%d}.A{ant:02d}.npz' if not dcbfile.is_file(): return with np.load(dcbfile,allow_pickle=True) as fid: dcb_gnss = fid['gnsses'] dcb_sol = fid['dcb_sol'] tecfiles = working_dir.glob(f'p05_*/{dt:%Y-%m-%d}/ROSW_TEC_{leo}_*_A{ant:02d}_*.nc') if len(out_tag) > 0: tar_dir = working_dir / 'p05_podTec' / f'{dt:%Y-%m-%d}_{out_tag}' tar_dir.mkdir(exist_ok=True,parents=True) with open('log/apply_dcb.txt','a') as fid: fid.write(f"{datetime.now():%H:%M:%S} - Start patching {iprof:4d}: ({dt:%Y-%m-%d}.{leo}.A{ant:02d})\n") for tecfile_orig in tecfiles: gnss = sub(r'.*_([A-Z]\d\d)_A.*',r'\1',tecfile_orig.name) if gnss not in dcb_gnss: continue dcb = dcb_sol[dcb_gnss==gnss][0] if out_tag == '': with Dataset(tecfile_orig, 'r+') as fid: fid['TEC'][...] = fid['TEC'][...] - fid.getncattr('leodcb') + dcb fid.setncattr('leodcb', dcb) fid.setncattr('tecmax', np.nanmax(fid['TEC'][...])) fid.setncattr('tecmin', np.nanmin(fid['TEC'][...])) else: with Dataset(tecfile_orig, 'r') as fid: dcb0 = fid.getncattr('leodcb') lv_err = fid.getncattr('leveling_err') tec_max = fid.getncattr('tecmax') -dcb0 + dcb tec_min = fid.getncattr('tecmin') -dcb0 + dcb tecfile = tar_dir / tecfile_orig.name if tecfile_orig == tecfile or lv_err>1 or tec_max>1000 or tec_min<-100: continue copy(tecfile_orig, tecfile) with Dataset(tecfile,'r+') as fid: fid['TEC'][...] = fid['TEC'][...] - fid.getncattr('leodcb') + dcb fid.setncattr('leodcb',dcb) fid.setncattr('tecmax', np.nanmax(fid['TEC'][...])) fid.setncattr('tecmin', np.nanmin(fid['TEC'][...])) with open('log/apply_dcb.txt','a') as fid: fid.write(f"{datetime.now():%H:%M:%S} - Done patching {iprof:4d}: ({dt:%Y-%m-%d}.{leo}.A{ant:02d})\n") return def process_apply_dcb(nProcess=10): dts = [datetime(2026, 1, 1) + timedelta(days=dd) for dd in range(263)] leos = ['C2E1', 'C2E2', 'C2E4', 'C2E5', 'C2E6'] ants = [1, 2] # for leo in leos: # for ant in ants: # dcbfile = f'p04_teczdcb/teczdcb_{leo}.{dts[-1]:%Y-%m-%d}.A{ant:02d}.npz' # with np.load(dcbfile,allow_pickle=True) as fid: # gnsses = fid['gnsses'] # dcb_sol = fid['dcb_sol'] # for dt in dts[-2::-1]: # dcbfile = f'p04_teczdcb/teczdcb_{leo}.{dt:%Y-%m-%d}.A{ant:02d}.npz' # with np.load(dcbfile,allow_pickle=True) as fid: # sec_tec = fid['sec_tec'] # gnsses0 = fid['gnsses'] # tecz = fid['tecz'] # dcb_sol0 = fid['dcb_sol'] # idx = ~np.isin(gnsses,gnsses0) # if np.any(idx): # dcb_sol1 = np.full_like(dcb_sol, np.nan) # dcb_sol1[~idx] = dcb_sol0 # dcb_sol0 = dcb_sol1 # daydiff = np.mean(dcb_sol[~idx] - dcb_sol0[~idx]) # dcb_sol0[idx] = dcb_sol[idx] - daydiff # np.savez(dcbfile, sec_tec=sec_tec, gnsses=gnsses, tecz=tecz, dcb_sol=dcb_sol0) # dcb_sol = dcb_sol0 allargs = np.array(np.meshgrid(dts, leos, ants, indexing='ij')).reshape((3,-1)).T allargs = list(np.hstack((np.arange(len(allargs)).reshape(-1,1),allargs))) import multiprocessing as mp with mp.Pool(processes=nProcess) as pool: pool.map(apply_dcb, allargs) return def calc_bore_sight(args): from re import sub import os from glob import glob from scipy.interpolate import CubicSpline from frame_conversions import ro_ecef2eci from traceback import extract_tb # c0 = 299792458 # kappa = 40.30819294518825 d_ref = args[0] leoName = args[1] ant_ref = args[2] gnss_ref = args[3] antennas = ['', 'POD1', 'POD2', 'OCC3', 'OCC4'] dt0_gps = datetime(1980, 1, 6) # ds_ref = [datetime.strptime(sub(r'.*/',r'',x),'%Y-%m-%d') for x in glob('./p01_podPair/*')] # ds_ref.sort() # # plt.close('all') # # fig1, axs1 = plt.subplots(3, 3, figsize=(12, 9)) # # fig2, axs2 = plt.subplots(3, 3, figsize=(12, 9)) (working_dir / 'p02_bsight').mkdir(exist_ok=True) out_dir = working_dir / f'./p02_bsight/{d_ref:%Y-%m-%d}/' out_dir.mkdir(exist_ok=True) # for leoName in leos_ref: # for gnss_ref in gnsses_ref: # for ant_ref in ants_ref: out_file = out_dir / f'bsight_{leoName}_{d_ref:%Y.%j}.A{ant_ref:02d}.{gnss_ref}.npz' if out_file.is_file(): with open("bore_sight_log.txt", 'a') as fid: fid.write(f"{datetime.now():%H:%M:%S} -- Skipping {out_file.stem}\n") return 1 reffiles = list(working_dir.glob(f'./p01_podPair/{d_ref:%Y-%m-%d}/occTab_{leoName}*A{ant_ref:02d}.{gnss_ref}.npz')) if len(reffiles) == 0: return -1 reffiles.sort() with open("bore_sight_log.txt", 'a') as fid: fid.write(f"{datetime.now():%H:%M:%S} -- Processing {out_file.stem}: {len(reffiles)} files.\n") mpe_0 = [] sec_ref_0 = [] sad_ref_0 = [] el_bsight_0 = [] az_bsight_0 = [] # print(f'{datetime.now():%H:%M:%S} - {d_ref:%Y-%m-%d} ({leoName}-{gnss_ref}-A{ant_ref:02d}), {len(reffiles):2d} files:', end='', flush=True) for reffile in reffiles: mpe = [] try: # print('.', end='', flush=True) dt_ref = [datetime.strptime(x, '%Y%m%d%H%M%S') for x in sub(r'.*(\d{14}\.\d{14}).*', r'\1', reffile.name).split('.')] dt_rr = datetime.strptime(sub(r'.*\.(\d{4}\.\d{3})\..*', r'\1', reffile.name), '%Y.%j') sec_rr = (dt_rr - dt0_gps) / timedelta(seconds=1) gnss_ref = sub(r'.*\.(...).npz', r'\1', reffile.name) ant_ref = int(sub(r'.*\.A(\d\d)\..*', r'\1', reffile.name)) with np.load('antenna_offset.npz') as fid: # antoffset_leo = fid[f'ANT_{leoName}_{antennas[ant_ref]}'][12:15] antbsight_leo = fid[f'ANT_{leoName}_{antennas[ant_ref]}'][15:18] # tecfiles = np.array(glob(f'/data3/xshao/planetiq/SWx/UCAR/L1B/20251220/podTec_GN05.*.{gnss_ref}.00*_nc')) # dts_ucar = np.array([datetime.strptime(sub(r'.*_GN0.\.(.*)\.\d{4}\.G.*',r'\1',tecfile),'%Y.%j.%H.%M') for tecfile in tecfiles]) # dur_ucar = np.array([timedelta(minutes=int(sub(r'.*\.(\d{4})\.G.*',r'\1',tecfile))) for tecfile in tecfiles]) # fidx = np.logical_and(dts_ucardt_ref[0]) # tecfiles = tecfiles[fidx] with np.load(reffile, allow_pickle=True) as fid: data_ref = fid['data'] header_ref = fid['header'] l_ref = np.where([x[0] == 'L' for x in header_ref])[0] s_ref = np.where([x[0] == 'S' for x in header_ref])[0] c_ref = np.where([x[0] == 'C' for x in header_ref])[0] sec_ref = (data_ref[:, 0] - sec_rr + data_ref[:, 1]).astype('int') f1_ref, f2_ref, _ = get_sat_freq(gnss_ref) try: Gdata, Gdt = read_gnss_v2(dt_ref, [gnss_ref]) except: continue tidx = np.logical_and(Gdt > dt_ref[0] - timedelta(minutes=10), Gdt < dt_ref[1] + timedelta(minutes=10)) dts_gnss = Gdt[tidx] sec_gnss = dt2float(dts_gnss, dt_rr).astype('int') pos_gnss = Gdata[tidx, 0, :4] vel_gnss = Gdata[tidx, 0, 4:] # OG = lagrange(sec_gnss/86400,Gdata[tidx, 0, :3],sec_ref/86400,11) # vel_gnss = Gdata[tidx, :, 4:] del Gdata, Gdt # Read in leo position/velocity/attitude files accordingly with np.load(f"p01_leoOrb/leoOrb_{leoName}.{dt_rr:%Y.%j}.npz", allow_pickle=True) as fid: Ldata = fid['data'] Ldts = fid['dts'] tidx = np.logical_and(Ldts > dt_ref[0] - timedelta(minutes=1), Ldts < dt_ref[1] + timedelta(minutes=1)) dts_leoorb = Ldts[tidx] pos_leo = Ldata[tidx, :4] vel_leo = Ldata[tidx, 4:] # sec_leo = dt2float(dts_leo, dt_rr).astype('int') # OL = lagrange(sec_leo/86400,Ldata[tidx, :3], sec_ref/86400, 11) # vel_leo = Ldata[tidx, 4:] del Ldata, Ldts with np.load(f"p01_leoAtt/leoAtt_{leoName}.{dt_rr:%Y.%j}.npz", allow_pickle=True) as fid: Ldata = fid['data'] Ldts = fid['dts'] Lvar = fid['vars'] tidx = np.logical_and(Ldts > dt_ref[0] - timedelta(minutes=1), Ldts < dt_ref[1] + timedelta(minutes=1)) if ~np.any(tidx): continue dts_leoatt = Ldts[tidx] sad_leo = Ldata[tidx, :][:, ['sad' in x for x in Lvar]] att_leo = Ldata[tidx, :][:, ['att' in x for x in Lvar]] aridx = np.cumsum(np.hstack(([False], np.sum(att_leo[1:] * att_leo[:-1], axis=1) < 0))) % 2 == 1 att_leo[aridx] *= -1 del Ldata, Ldts, Lvar tidx = np.logical_and(dts_leoorb >= dts_leoatt[0], dts_leoorb <= dts_leoatt[-1]) if ~np.any(tidx): continue dts_leo = dts_leoorb[tidx] sec_leo = dt2float(dts_leo, dt_rr) pos_leo = pos_leo[tidx] vel_leo = vel_leo[tidx] dts_leoatt, tidx = np.unique(dts_leoatt, return_index=True) # sad_leo = CubicSpline(dt2float(dts_leoatt, dts_leo[0]), sad_leo[tidx])(dt2float(dts_leo, dts_leo[0])) # att_leo = CubicSpline(dt2float(dts_leoatt, dts_leo[0]), att_leo[tidx])(dt2float(dts_leo, dts_leo[0])) att_leo /= np.sum(att_leo ** 2, axis=1)[:, None] ** 0.5 # del dts_leoorb, dts_leoatt if sec_leo[0] > sec_ref[0] or sec_leo[-1] < sec_ref[-1]: # raise ValueError( # f'LEO time ({sec_leo[0]:.1f} >> {sec_leo[-1]:.1f}) could not cover REF event ({sec_ref[0]:.1f} >> {sec_ref[-1]:.1f})!') continue pos_leo_eci, vel_leo_eci = ro_ecef2eci(pos_leo, vel_leo, dts_leo) pos_ref_eci, vel_ref_eci = ro_ecef2eci(pos_gnss[:, :3], vel_gnss, dts_gnss) OL = lagrange(sec_leo / 86400, pos_leo_eci, sec_ref / 86400, 11) OG = lagrange(sec_gnss / 86400, pos_ref_eci, sec_ref / 86400, 11) antbsight_eci = inv_quat_rotate(att_leo, antbsight_leo) ex = lagrange(dt2float(dts_leoatt, dt_rr) / 86400, antbsight_eci, sec_ref / 86400, 11) ex /= np.sum(ex ** 2, axis=1)[:, None] ** 0.5 ez = -OL # * [1, 1, 1.00669437999014] ey = np.cross(ez, ex) ey /= np.sum(ey ** 2, axis=1)[:, None] ** 0.5 ez = np.cross(ex, ey) LG = OG - OL LG_bsight = np.array([np.sum(LG * x, axis=1) for x in [ex, ey, ez]]).T el_bsight = -np.arcsin(LG_bsight[:, 2] / np.sum(LG_bsight ** 2, axis=1) ** 0.5) / np.pi * 180 az_bsight = np.arctan2(LG_bsight[:, 1], LG_bsight[:, 0]) / np.pi * 180 sad_leo -= np.cumsum(np.round(np.diff(sad_leo, axis=0, prepend=sad_leo[:1]) / 360) * 360, axis=0) sad_ref = lagrange(dt2float(dts_leoatt, dt_rr) / 86400, sad_leo, sec_ref / 86400, 1) sad_ref[:,0] = sad_ref[:,0] % 360 sad_ref[:,1] = (sad_ref[:,1] + 180) % 360 - 180 # idx_sad = np.isin(np.round(sad_ref * 10), np.arange(-1200, 1500, 300) - 123) proj_lg = OL - LG * (np.sum(OL * LG, axis=1) / np.sum(LG * LG, axis=1))[:, None] # altp_km = EarthLocation.from_geocentric(*proj_lg.T, unit='km').geodetic.height.value # altl_km = EarthLocation.from_geocentric(*OL.T, unit='km').geodetic.height.value # coszn = np.sum(OL * LG, axis=1) / np.sum(LG * LG, axis=1) ** 0.5 / np.sum(OL * OL, axis=1) ** 0.5 # Pdiff = data_ref[:,l_ref[0]] - data_ref[:,l_ref[1]] fr_ref = 2 / ((f1_ref / f2_ref) ** 2 - 1) + 1 S = data_ref[:, s_ref] S2 = S ** 2 # cumS = np.cumsum(np.vstack(([[0, 0]], S)), axis=0) # cumS2 = np.cumsum(np.vstack(([[0, 0]], S2)), axis=0) # clen = 31 # S4 = ((cumS2[clen:] - cumS2[:-clen]) * clen / (cumS[clen:] - cumS[:-clen]) ** 2 - 1) ** .5 # L = data_ref[:, l_ref] - data_ref[:, c_ref] # L2 = L ** 2 # cumL = np.cumsum(np.vstack(([[0, 0]], L)), axis=0) # cumL2 = np.cumsum(np.vstack(([[0, 0]], L2)), axis=0) # sigma = ((cumL2[clen:] - cumL2[:-clen]) * clen / (cumL[clen:] - cumL[:-clen]) ** 2 - 1) ** .5 S2 /= np.sum(S2, axis=0)[None,:] D1 = data_ref[:, c_ref[0]] - fr_ref * data_ref[:, l_ref[0]] + (fr_ref - 1) * data_ref[:, l_ref[1]] D2 = data_ref[:, c_ref[1]] - (fr_ref + 1) * data_ref[:, l_ref[0]] + fr_ref * data_ref[:, l_ref[1]] sig1 = np.sum(D1 * S2[:, 0]) sig2 = np.sum(D2 * S2[:, 1]) # sad_ref = np.round(sad_ref * 10) # idx_sad = np.isin(sad_ref, np.arange(-1200, 1500, 300) - 123) # sad_ref /= 10 idx_sad = np.logical_and(sad_ref[:,1]>-180, sad_ref[:,1]<180) mpe = np.array([D1[idx_sad] - sig1, D2[idx_sad] - sig2]).T sec_ref = sec_ref[idx_sad] sad_ref = sad_ref[idx_sad] el_bsight = el_bsight[idx_sad] az_bsight = az_bsight[idx_sad] except Exception as e: with open("bore_sight_log.txt", 'a') as fid: fid.write(f"{datetime.now().strftime('%H:%M:%S')} -- Error detected for {reffile}:\n") fid.write(f' --- {type(e).__name__} >> {str(e)}\n') for tb in extract_tb(e.__traceback__): fid.write(f' --- line {tb.lineno} of {tb.filename}\n') if len(mpe) == 0: continue mpe_0.append(mpe) sec_ref_0.append(sec_ref) sad_ref_0.append(sad_ref) el_bsight_0.append(el_bsight) az_bsight_0.append(az_bsight) if len(mpe_0) > 1: mpe_0 = np.vstack(mpe_0) sec_ref_0 = np.hstack(sec_ref_0) sad_ref_0 = np.vstack(sad_ref_0) el_bsight_0 = np.hstack(el_bsight_0) az_bsight_0 = np.hstack(az_bsight_0) elif len(mpe_0) == 1: mpe_0 = mpe_0[0] sec_ref_0 = sec_ref_0[0] sad_ref_0 = sad_ref_0[0] el_bsight_0 = el_bsight_0[0] az_bsight_0 = az_bsight_0[0] else: # print('', flush=True) return -1 sec_ref_0, idx_sec = np.unique(sec_ref_0, return_index=True) mpe_0 = mpe_0[idx_sec] sad_ref_0 = sad_ref_0[idx_sec] el_bsight_0 = el_bsight_0[idx_sec] az_bsight_0 = az_bsight_0[idx_sec] np.savez(out_file, mpe=mpe_0, sec_ref=sec_ref_0, sad_ref=sad_ref_0, el_bsight=el_bsight_0, az_bsight=az_bsight_0) # print('', flush=True) return 0 def process_calc_bore_sight(nProcess=20): leos = ['C2E1','C2E2','C2E4','C2E5','C2E6'] ants = [1,2] gnsses = [f'G{ii+1:02d}' for ii in range(32)]+[f'R{ii+1:02d}' for ii in range(28)] ds_ref = [datetime(2025,1,1) + timedelta(days=x) for x in range(575)] allargs = [(d_ref, leo, ant, gnss) for d_ref in ds_ref for leo in leos for ant in ants for gnss in gnsses] if nProcess > 0: import multiprocessing as mp with mp.Pool(processes=nProcess) as pool: pool.map(calc_bore_sight, allargs) else: for args in allargs: calc_bore_sight(args) def gen_bore_sight_lut(): from glob import glob from re import sub import pickle degsize = 2 el_rng = [-90, 90] az_rng = [-180, 180] el_size = int((el_rng[1] - el_rng[0]) // degsize) az_size = int((az_rng[1] - az_rng[0]) // degsize) el_bins = np.linspace(el_rng[0], el_rng[1], el_size + 1) az_bins = np.linspace(az_rng[0], az_rng[1], az_size + 1) with open('p01_podCode/podCode_Psats.pkl', 'rb') as fid: ptimes = pickle.load(fid) # for leoName in ['C2E1', 'C2E2', 'C2E4', 'C2E5', 'C2E6']: for leoName in ['C2E2']: ll = int(leoName[3]) - 1 for ant_ref in [2]: for gnss_ref in ['GL1','GC2','GP2','R01','R02','R03','R04','R09','R10','R11','R12','R17','R18','R19','R20']: bsfiles = glob(f'./p02_bsight/*/bsight_{leoName}*.A{ant_ref:02d}.{gnss_ref[0]}*.npz') bsfiles.sort() bsgnsses = [sub(r'.*([A-Z]\d\d)\.npz', r'\1', bsfile) for bsfile in bsfiles] bsdates = [datetime.strptime(sub(r'.*_(\d{4}\.\d{3})\..*', r'\1', bsfile), '%Y.%j') for bsfile in bsfiles] # bsfiles = [x for x in bsfiles if datetime.strptime(sub(r'.*_(\d{4}\.\d{3})\..*', r'\1', x), '%Y.%j') <= datetime(2026,1,20)] date_rng = [datetime.strptime(sub(r'.*_(\d{4}\.\d{3})\..*', r'\1', bsfiles[x]), '%Y.%j') for x in [0,-1]] gnss_fname = gnss_ref outfile = f'./p03_bslut/bslut_{leoName}.{date_rng[0]:%Y-%m-%d}.D{int((date_rng[1]-date_rng[0])/timedelta(days=1)+1):03d}.A{ant_ref:02d}.{gnss_fname}.npz' print(f'{datetime.now():%H:%M:%S} - {outfile}', flush=True) if Path(outfile).is_file(): continue mpe = [] eat_dscr = [] sda_dscr = [] for bsfile, bsgnss, bsdate in zip(bsfiles, bsgnsses, bsdates): if gnss_ref[0] == 'G': if gnss_ref == 'GC2': if bsgnss in ptimes and (ptimes[bsgnss][ll] == None or ptimes[bsgnss][ll] > bsdate + timedelta(days=1)): continue if gnss_ref == 'GP2': if bsgnss not in ptimes or (ptimes[bsgnss][ll] != None and ptimes[bsgnss][ll] < bsdate): continue elif gnss_ref[0] == 'R': bsgnum = int(bsgnss[1:]) rfgnum = int(gnss_ref[1:]) if bsgnss[0] != 'R' or bsgnum % 4 != rfgnum % 4 or (bsgnum-1) // 8 != (rfgnum-1) // 8: continue with np.load(bsfile, allow_pickle=True) as fid: mpe.append(fid['mpe'][...]) sad_ref0 = fid['sad_ref'][...] el_bsight0 = fid['el_bsight'][...] az_bsight0 = fid['az_bsight'][...] idx = np.all((sad_ref0[:, 0] >= 0, sad_ref0[:, 0] <= 360, sad_ref0[:, 1] >= -45, sad_ref0[:, 1] < 45), axis=0) el_dscr = ((el_bsight0[idx] - el_rng[0]) // degsize).astype('int') el_dscr[el_dscr == el_size] -= 1 az_dscr = (((az_bsight0[idx] - az_rng[0]) % (az_rng[1] - az_rng[0])) // degsize).astype('int') sta_dscr = ((sad_ref0[idx, 1] + 45) // 30).astype('int') eat_dscr.append((el_dscr * az_size + az_dscr) * 3 + sta_dscr) sda_dscr.append((((sad_ref0[idx, 0] + 185) % 360) // 10).astype('int')) mpe = np.vstack(mpe) eat_dscr = np.hstack(eat_dscr) sda_dscr = np.hstack(sda_dscr) # idx = np.all(np.abs(mpe) < 3, axis=1) # mpe = mpe[idx] # sad_ref = sad_ref[idx] # el_bsight = el_bsight[idx] # az_bsight = az_bsight[idx] sda_bins = np.unique(sda_dscr) mpe_grid = np.full((el_size, az_size, 3, 36, 2), np.nan) idx = np.lexsort((sda_dscr,eat_dscr)) mpe = mpe[idx] eat_dscr = eat_dscr[idx] sda_dscr = sda_dscr[idx] div = np.hstack(([0],np.where(np.diff(eat_dscr)>0)[0]+1,[len(eat_dscr)])) for dd in range(len(div)-1): if dd%1000==0: print(dd) el_idx = int(eat_dscr[div[dd]] // 3 // az_size) az_idx = int(eat_dscr[div[dd]] // 3) % az_size sta_idx = int(eat_dscr[div[dd]] % 3) sda_slice = sda_dscr[div[dd]:div[dd+1]] mpe_slice = mpe[div[dd]:div[dd+1]] sdiv = np.hstack(([0],np.where(np.diff(sda_slice)>0)[0]+1,[len(sda_slice)])) for ss in range(len(sdiv)-1): if sdiv[ss + 1] - sdiv[ss] < 3: continue sda_idx = sda_slice[sdiv[ss]] for bb in range(2): mpe_bin = mpe_slice[sdiv[ss]:sdiv[ss + 1], bb] mpe_bin = mpe_bin[np.abs(mpe_bin)<3] mpe_med = np.median(mpe_bin) mpe_std = np.std(mpe_bin) while np.any(np.abs(mpe_bin - mpe_med) > 3 * mpe_std): mpe_bin = mpe_bin[np.abs(mpe_bin - mpe_med) < 3 * mpe_std] mpe_med = np.median(mpe_bin) mpe_std = np.std(mpe_bin) mpe_grid[el_idx, az_idx, sta_idx, sda_idx, bb] = np.mean(mpe_bin) for bb in range(2): sda_idx = np.where(np.isfinite(mpe_grid[el_idx, az_idx, sta_idx, :, bb]))[0] if len(sda_idx) > 1: mpe_grid[el_idx, az_idx, sta_idx, :, bb] = np.interp(np.arange(36), sda_idx, mpe_grid[el_idx, az_idx, sta_idx, sda_idx, bb], left=np.nan, right=np.nan) np.savez(outfile, el_bins=el_bins, az_bins=az_bins, sad_bins=sda_bins, mpe_grid=mpe_grid) return def patch_edp(dt): from subprocess import run from re import sub print(f'{datetime.now():%H:%M:%S} - Processing {dt:%Y-%m-%d}...',flush=True) src_dir = working_dir / f'p05_podTec/{dt:%Y-%m-%d}' tar_dir = working_dir / f'p06_podEdp/{dt:%Y-%m-%d}' tar_dir.mkdir(parents=True, exist_ok=True) for cnt, src_file in enumerate(src_dir.glob('*.nc')): # if cnt%50==0: # print(f'{datetime.now():%H:%M:%S} - ',end='',flush=True) tar_file = tar_dir / sub(r'_TEC_',r'_EDP_',src_file.name) if tar_file.is_file(): # print('-', end='',flush=True) continue results = run([(working_dir / '.venv/bin/python').as_posix(), 'TEC_to_EDP_final_gnss.py', '--TEC', src_file.as_posix(), '--out', tar_dir.as_posix()], capture_output=True, cwd=working_dir.parent / 'TEC_EDP') # if tar_file.is_file(): # print('o', end='',flush=True) # else: # print('x', end='',flush=True) # if cnt%50==49: # print(flush=True) # print(flush=True) print(f'{datetime.now():%H:%M:%S} - Done {dt:%Y-%m-%d}...', flush=True) return