import numpy as np from re import sub # from glob import glob from astropy.coordinates import EarthLocation from traceback import extract_tb from datetime import datetime, timedelta from pathlib import Path from tools import dt2float from netCDF4 import Dataset working_dir = Path(__file__).absolute().parent def log_message(msg, working_leo=''): with open(working_dir / f"log/progress_dcb_{working_leo}_{datetime.now():%Y%m}.txt", 'a') as fid: fid.write(f"{datetime.now():%H:%M:%S} {msg}\n") alti_km = 2400 dt0_gps = datetime(1980,1,6) leos = ['GN04', 'GN05'] # TC: updated 2026/05/18 (1 line, fix sign typo) dts = [datetime.today() + timedelta(days=dd) for dd in range(-1,0)] # dts = [datetime(2026,1,1) + timedelta(days=dd) for dd in range(1)] ants = [0, 1] for dt in dts: sec_rr = (dt - dt0_gps) / timedelta(seconds=1) for ll, leo in enumerate(leos): for aa, ant in enumerate(ants): log_message(f"- Processing ({dt:%Y-%m-%d}.{leo}.A{ant:02d})") try: outfile = working_dir / f'p04_teczdcb/teczdcb_{leo}.{dt:%Y-%m-%d}.A{ant:02d}.npz' for tdif in range(1,366): oldfile = working_dir / f'p04_teczdcb/teczdcb_{leo}.{dt-timedelta(days=tdif):%Y-%m-%d}.A{ant:02d}.npz' if oldfile.exists(): with np.load(oldfile, allow_pickle=True) as fid: gnsses = fid['gnsses'] dcb_sol = fid['dcb_sol'] break alltec = np.full((187200, len(gnsses)), np.nan) allcoszn = np.full((187200, len(gnsses)), np.nan) allmeps = np.full((187200, len(gnsses)), np.nan) orbfile = working_dir / f'p01_leoOrb/leoOrb_{leo}.{dt-timedelta(days=1):%Y.%j}.npz' with np.load(orbfile, allow_pickle=True) as fid: Ldata0 = fid['data'] Ldts0 = fid['dts'] orbfile = working_dir / f'p01_leoOrb/leoOrb_{leo}.{dt:%Y.%j}.npz' # TC: updated 2026/05/18 (10 lines, check existence of next day orbit file) if orbfile.exists(): with np.load(orbfile, allow_pickle=True) as fid: Ldata1 = fid['data'] Ldts1 = fid['dts'] idx = ~np.isin(dt2float(Ldts0,dt),dt2float(Ldts1,dt),assume_unique=True) Ldts = np.hstack((Ldts0[idx], Ldts1)) Ldata = np.vstack((Ldata0[idx], Ldata1)) else: Ldts = Ldts0 Ldata = Ldata0 idx = np.argsort(Ldts) Ldts = Ldts[idx] Ldata = Ldata[idx] sec_ref = np.round(((Ldts - dt) / timedelta(seconds=1)).astype('float')).astype('int') idx = sec_ref>=-93600 leo_pos = EarthLocation(x=Ldata[idx,0],y=Ldata[idx,1],z=Ldata[idx,2],unit='km') alllon, alllat, _ = leo_pos.to_geodetic() allhr = (sec_ref[idx]/3600+alllon.value/15+12)%24-12 tecfiles = list( working_dir.glob(f'p05_*/{dt - timedelta(days=1):%Y-%m-%d}/ROSW_TEC_{leo}_*_A{ant:02d}_*.nc')) + list( working_dir.glob(f'p05_*/{dt:%Y-%m-%d}/ROSW_TEC_{leo}_*_A{ant:02d}_*.nc')) for tecfile in tecfiles: try: gnss = sub(r'.*_([A-Z]\d\d)_A\d\d_.*',r'\1',tecfile.name) if gnss not in gnsses: gnsses = np.append(gnsses, gnss) dcb_sol = np.append(dcb_sol, np.nan) alltec = np.hstack((alltec, np.full_like(alltec[:,:1], np.nan))) allcoszn = np.hstack((allcoszn, np.full_like(allcoszn[:, :1], np.nan))) allmeps = np.hstack((allmeps, np.full_like(allmeps[:, :1], np.nan))) gidx = list(gnsses).index(gnss) with Dataset(tecfile, 'r') as fid: tidx = (fid['time'][...]-sec_rr).astype('int')+93600 dcb = fid.getncattr('leodcb') tec = fid['TEC'][...]-dcb coszn = np.sin(fid['elevation'][...]*np.pi/180) OL = np.column_stack([fid['x_LEO'][...],fid['y_LEO'][...],fid['z_LEO'][...]]) altl_km = EarthLocation.from_geocentric(*OL.T, unit='km').geodetic.height.value rad_ratio = (alti_km + 6371) / (altl_km + 6371) meps = (coszn + (rad_ratio ** 2 + coszn ** 2 - 1) ** 0.5) / (1 + rad_ratio) alltec[tidx, gidx] = tec allcoszn[tidx, gidx] = coszn allmeps[tidx, gidx] = meps except: continue mask = allcoszn > .5 # tmask = np.full((187200,), False) # tmask[sec_ref[idx]+7200] = np.logical_and(np.abs(allhr) < 4, np.abs(alllat.value) > 30) # nprof = np.sum(mask, axis=1) T = alltec.copy() T[~mask] = np.nan M = allmeps.copy() M[~mask] = np.nan gidx = np.where(np.any(np.isfinite(T), axis=0))[0] # allcoszn = allcoszn[:, gidx] # alltec = alltec[:, gidx] # allmeps = allmeps[:, gidx] # M = M[:, gidx] # T = T[:, gidx] # mask = allcoszn > .5 # allcoszn0 = allcoszn.copy() # alltec0 = alltec.copy() # allmeps0 = allmeps.copy() # gidx0 = gidx.copy() dcb_sol_0 = dcb_sol.copy() modtec = (T + dcb_sol_0) * M idx = np.any((modtec > 50, modtec < 0, np.isnan(modtec)), axis=0) T[idx] = np.nan M[idx] = np.nan T0 = T.copy() M0 = M.copy() idx = np.full((T.shape[1],),False) while True: dcb_sol[:] = np.nan T = T0.copy() M = M0.copy() T[:,idx] = np.nan M[:,idx] = np.nan modtec = (T + dcb_sol_0) * M idx = np.abs(modtec - np.nanmedian(modtec, axis=1)[:, None]) > 5 T[idx] = np.nan M[idx] = np.nan pidx = np.logical_and(np.isfinite(T),np.isfinite(M)) A = np.zeros((T.shape[1], T.shape[1])) B = np.zeros((T.shape[1], )) for t in range(T.shape[0]): idx = pidx[t] nprof = np.sum(idx) if nprof < 2: continue A[idx,idx] += M[t,idx]**2 A[np.ix_(idx, idx)] -= M[t,idx,None]*M[None,t,idx]/nprof B[idx] += (M[t,idx]@T[t,idx])*M[t,idx]/nprof - T[t,idx]*M[t,idx]**2 idx = np.logical_or(np.all(A==0,axis=0),np.all(A==0,axis=1)) if np.any(idx): A = A[~idx][:,~idx] B = B[~idx] gidx = np.where(~idx)[0] # gidx = gidx[~idx] # allcoszn = allcoszn[:, ~idx] # alltec = alltec[:, ~idx] # allmeps = allmeps[:, ~idx] # M = M[:, ~idx] # T = T[:, ~idx] # mask = allcoszn > .5 if np.linalg.matrix_rank(A) 3 while np.any(idx): idx = np.abs(modtec - np.nanmedian(modtec, axis=1)[:, None]) > 2.9 dcb_sol[:] = np.nan T[idx] = np.nan M[idx] = np.nan pidx =np.logical_and(np.isfinite(T),np.isfinite(M)) A = np.zeros((T.shape[1], T.shape[1])) B = np.zeros((T.shape[1],)) for t in range(T.shape[0]): idx = pidx[t] nprof = np.sum(idx) if nprof < 2: continue A[idx, idx] += M[t, idx] ** 2 A[np.ix_(idx, idx)] -= M[t, idx, None] * M[None, t, idx] / nprof B[idx] += (M[t, idx] @ T[t, idx]) * M[t, idx] / nprof - T[t, idx] * M[t, idx] ** 2 idx = np.logical_or(np.all(A==0,axis=0),np.all(A==0,axis=1)) if np.any(idx): A = A[~idx][:, ~idx] B = B[~idx] gidx = np.where(~idx)[0] # gidx = gidx[~idx] # allcoszn = allcoszn[:, ~idx] # alltec = alltec[:, ~idx] # allmeps = allmeps[:, ~idx] # M = M[:, ~idx] # T = T[:, ~idx] # mask = allcoszn > .1 if np.linalg.matrix_rank(A) < A.shape[0]: raise ValueError(f"Singular linear system (2) ({dt:%Y-%m-%d}.{leo}.A{ant:02d})") dcb_sol[gidx] = np.linalg.solve(A, B) modtec = (T + dcb_sol) * M # modtec[~mask] = np.nan idx = np.abs(modtec - np.nanmedian(modtec, axis=1)[:, None]) > 3 idx = np.logical_and(np.abs(dcb_sol - dcb_sol_0) > (100 if dt == datetime(2026, 1, 1) else 1), np.abs(dcb_sol - dcb_sol_0 - np.nanmedian(dcb_sol - dcb_sol_0)) > 0.75) if ~np.any(idx): break else: log_message(f'-- bad solution for ({dt:%Y-%m-%d}.{leo}.A{ant:02d})') for ii in np.where(idx)[0]: log_message(f'--- {gnsses[ii]}:{dcb_sol_0[ii]:8.2f} -> {dcb_sol[ii]:8.2f}') idx = np.logical_or(idx,np.isnan(dcb_sol)) # idx = np.all((np.isfinite(dcb_sol),np.isfinite(dcb_sol_0),np.abs(dcb_sol-dcb_sol_0)>(100 if dt==dts[0] else 1)),axis=0) # if np.any(idx): # dcb_sol[idx] = np.nan # modtec = (T + dcb_sol) * M # modtec[~mask] = np.nan tecz = np.nanmean(modtec, axis=1) idx = np.logical_and(np.isfinite(dcb_sol),np.isfinite(dcb_sol_0)) daydiff = 0 if np.any(idx): daydiff = np.mean(dcb_sol[idx]-dcb_sol_0[idx]) for gg in range(len(dcb_sol)): if np.isfinite(dcb_sol[gg]): continue dcb0 = dcb_sol_0[gg] + daydiff if np.isfinite(dcb0): T = T0[:,gg] M = M0[:,gg] idx = np.all(np.isfinite([T,M,tecz]),axis=0) dcb = tecz[idx]/M[idx]-T[idx] idx = np.abs((T[idx]+dcb0)*M[idx]-tecz[idx])<3 else: idx = allcoszn[:,gg]>0.5 T = alltec[:,gg] M = allmeps[:,gg] T[~idx] = np.nan M[~idx] = np.nan idx = np.all(np.isfinite([T, M, tecz]), axis=0) dcb = tecz[idx] / M[idx] - T[idx] idx = np.abs((T[idx]+np.median(dcb))*M[idx]-tecz[idx])<3 if np.sum(idx)>6000: dcb_sol[gg] = np.nanmean(dcb[idx]) if np.isfinite(dcb0): diff = min(max(dcb_sol[gg]-dcb0,-0.2),0.2) dcb_sol[gg] = dcb0+diff elif np.isfinite(dcb0): dcb_sol[gg] = dcb0 gidx = np.isfinite(dcb_sol) dcb_sol = dcb_sol[gidx] gnsses = gnsses[gidx] gidx = np.argsort(gnsses) dcb_sol = dcb_sol[gidx] gnsses = gnsses[gidx] np.savez(outfile, sec_tec=np.arange(-93600,93600),gnsses=gnsses,tecz=tecz,dcb_sol=dcb_sol) # sTM2 = np.nansum(T*M**2, axis=1) # sM = np.nansum(M, axis=1) # sM2 = np.nansum(M**2, axis=1) # sTM = np.nansum(T*M, axis=1) # A = sM2-sM**2/nprof # B = -sTM2+sM*sTM/nprof # A[nprof<2] = np.nan # B[nprof<2] = np.nan # leodcb = (B/A)[:,None] # modtec = (alltec+leodcb)*allmeps # modtec[~mask] = np.nan # idx = np.logical_and(np.nanstd(modtec,axis=1)<0.1,np.nanmean(modtec,axis=1)>0) log_message(f"- Done ({dt:%Y-%m-%d}.{leo}.A{ant:02d})") except Exception as e: log_message(f'-- Error detected for ({dt:%Y-%m-%d}.{leo}.A{ant:02d})') log_message(f'--- {type(e).__name__} >> {str(e)}') for tb in extract_tb(e.__traceback__): log_message(f"--- line {tb.lineno} of {tb.filename}")