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 if '__file__' in dir() else Path.cwd() 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") def process_dcb(args): leo = args[0] ant = args[1] dt0 = args[2] dur = args[3] alti_km = 2400 dt0_gps = datetime(1980,1,6) # TC: updated 2026/05/18 (1 line, fix sign typo) # dts = [datetime.today() + timedelta(days=dd) for dd in range(-1,0)] dts = [dt0 + timedelta(days=dd) for dd in range(dur)] for dt in dts: sec_rr = (dt - dt0_gps) / timedelta(seconds=1) 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.is_file(): with np.load(oldfile, allow_pickle=True) as fid: gnsses = fid['gnsses'] dcb_sol = fid['dcb_sol'] break if not oldfile.is_file(): gnsses = np.array([]) dcb_sol = np.array([]) 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((alltec.shape[0],1), np.nan))) allcoszn = np.hstack((allcoszn, np.full((allcoszn.shape[0],1), np.nan))) allmeps = np.hstack((allmeps, np.full((allmeps.shape[0],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'][...]]) lv_err = fid.getncattr('leveling_err') tec_max = fid.getncattr('tecmax') tec_min = fid.getncattr('tecmin') if lv_err>1 or tec_max>1500 or tec_min<-500: continue 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 = np.logical_and(allcoszn > np.cos(80/180*np.pi), allcoszn < np.cos(40/180*np.pi)) # 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() if np.any(np.isfinite(dcb_sol_0)): modtec = (T + dcb_sol_0) * M idx = np.any((modtec > 100, modtec < -50, np.isnan(modtec)), axis=0) T[idx] = np.nan M[idx] = np.nan T0 = T.copy() M0 = M.copy() # gidx = np.array([gnss[0]=='G' for gnss in gnsses]) # T0[:,~gidx] = np.nan # M0[:,~gidx] = np.nan 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 if np.any(np.isfinite(dcb_sol_0)): modtec = (T + dcb_sol_0) * M idx = np.abs(modtec - np.nanmedian(modtec, axis=1)[:, None]) > 3 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 weight = 1 A[idx, idx] += M[t, idx] ** 2 * weight A[np.ix_(idx, idx)] -= M[t, idx, None] * M[None, t, idx] / nprof * weight B[idx] += ((M[t, idx] @ T[t, idx]) * M[t, idx] / nprof - T[t, idx] * M[t, idx] ** 2) * weight 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 A.shape[0]==0 or 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 weight = 1 A[idx, idx] += M[t, idx] ** 2 * weight A[np.ix_(idx, idx)] -= M[t, idx, None] * M[None, t, idx] / nprof * weight B[idx] += ((M[t, idx] @ T[t, idx]) * M[t, idx] / nprof - T[t, idx] * M[t, idx] ** 2) * weight 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 A.shape[0]==0 or 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 5), np.abs(dcb_sol - dcb_sol_0 - np.nanmedian(dcb_sol - dcb_sol_0)) > 5) 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 # T0 = alltec.copy() # M0 = allmeps.copy() # T_shift = T0 + dcb_sol # idx = np.logical_and(np.abs(np.diff(T_shift, axis=0, prepend=T_shift[:1])) < 1, # np.abs(np.diff(T_shift, axis=0, append=T_shift[-1:])) < 1) # T_shift[~idx] = np.nan # dcb_shift = -np.nanmin(T_shift) # dcb_sol += dcb_shift # print(f'{leo}_A{ant:02d}: {dcb_shift}') # modtec = (T + dcb_sol) * M 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] modtec = (T[idx] + dcb0) * M[idx] idx = np.logical_and(np.abs(modtec - tecz[idx]) < 3, modtec >= 0) if not np.isfinite(dcb0) or np.sum(idx) < 1000: idx = np.logical_and(allcoszn[:,gg] > np.cos(80/180*np.pi), allcoszn[:,gg] < np.cos(40/180*np.pi)) 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] modtec = (T[idx] + np.median(dcb)) * M[idx] idx = np.logical_and(np.abs(modtec - tecz[idx]) < 3, modtec >= 0) if np.sum(idx)>6000: dcb_sol[gg] = np.nanmean(dcb[idx]) if np.isfinite(dcb0): driftbound = 0.5 if gnsses[gg].startswith('C') else 0.2 diff = min(max(dcb_sol[gg]-dcb0,-driftbound),driftbound) 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}") if __name__ == '__main__': leos = ['C2E1', 'C2E2', 'C2E4', 'C2E5', 'C2E6'] ants = [1, 2] leos, ants = np.meshgrid(leos, ants, indexing='ij') leos = leos.ravel() ants = ants.ravel() # dt0 = [datetime(2025,12,31)] * len(leos) # dur = [264] * len(leos) dt0 = [datetime.now().replace(hour=0, minute=0, second=0, microsecond=0) - timedelta(days=1)] * len(leos) dur = [1] * len(leos) allargs = list(zip(leos,ants,dt0,dur)) nProcess = len(allargs) import multiprocessing as mp with mp.Pool(processes=nProcess) as pool: pool.map(process_dcb, allargs)