import numpy as np from datetime import datetime, timedelta # from glob import glob # import os from pathlib import Path from tools import to_datetime from re import sub working_dir = Path(__file__).absolute().parent lsat_names = dict(GN04=4, GN05=5, YM08=8) for keys in list(lsat_names.keys()): lsat_names[lsat_names[keys]] = keys def log_message(msg, working_leo): with open(Path(__file__).parent / f"log/progress_rt_{working_leo}_{datetime.now():%Y%m%d}.txt", 'a') as fid: fid.write(f"{datetime.now():%H:%M:%S} {msg}\n") def extract_att(working_leo): att_dir = working_dir / f'L1_extract_{working_leo}/leoAtt' attfiles = list(att_dir.glob('*_txt')) leo_sats = np.unique([lsat_names[int(sub(r'.*_\d{4}\.\d{3}\.(\d{3})\.\d\d_.*', r'\1', x.name))] for x in attfiles]) # attfiles = glob(f'/data3/xshao/planetiq/SWx/L1A/*/gnomes/level1a/leoAtt/leoAtt_*') # leo_sats = np.unique([lsat_names[int(sub(r'.*_\d{4}\.\d{3}\.(\d{3})\.\d\d_.*', r'\1', x))] for x in attfiles]) # leo_dts = np.unique([datetime.strptime(sub(r'.*_(\d{4}\.\d{3})\..*',r'\1',x),'%Y.%j') for x in attfiles]) vars_list = ['att_x', 'att_y', 'att_z', 'att_w', 'ang_x', 'ang_y', 'ang_z', 'pe_x', 'pe_y', 'pe_z', 've_x', 've_y', 've_z', 'pi_x', 'pi_y', 'pi_z', 'vi_x', 'vi_y', 'vi_z', 'thf', 'sad_a','sad_b1','sad_b2'] for ff in leo_sats: attfiles = list(att_dir.glob(f'*.*.{lsat_names[ff]:03d}.*_txt')) attfiles.sort() dt = datetime.strptime(sub(r'.*_(\d{4}\.\d{3})\.\d{3}\.\d\d_.*', r'\1', attfiles[0].name), '%Y.%j') attfileout = working_dir / f'p01_leoAtt/leoAtt_{ff}.{dt:%Y.%j}.npz' attfileprev = working_dir / f'p01_leoAtt/leoAtt_{ff}.{dt-timedelta(days=1):%Y.%j}.npz' dts0 = [] data0 = [] tmp0 = np.full((len(vars_list),),np.nan) for attfile in attfiles: for line in open(attfile,'r'): if line.startswith("tim"): if len(dts0)>len(data0): data0.append(tmp0.copy()) tmp0[:] = np.nan dts0.append(to_datetime(line[3:].split())) elif line.startswith("att"): tmp0[:4] = np.fromstring(line[3:], sep=" ")[1:-1] elif line.startswith("ang"): tmp0[4:7] = np.fromstring(line[3:], sep=" ") elif line.startswith("pve"): tmp0[7:13] = np.fromstring(line[3:], sep=" ") elif line.startswith("pvi"): tmp0[13:19] = np.fromstring(line[3:], sep=" ") elif line.startswith("thf"): tmp0[19] = np.fromstring(line[3:], sep=" ") elif line.startswith("sad"): tmp = np.fromstring(line[3:], sep=" ") tmp0[20:20+len(tmp)] = tmp data0.append(tmp0.copy()) dts0 = np.array(dts0) data0 = np.array(data0) if attfileout.exists(): with np.load(attfileout, allow_pickle=True) as fid: data = fid['data'] dts = fid['dts'] idx = ~np.isin(dts, dts0) dts = np.hstack([dts[idx],dts0]) data = np.vstack([data[idx], data0]) idx = np.argsort(dts) dts = dts[idx] data = data[idx] elif attfileprev.exists(): with np.load(attfileprev, allow_pickle=True) as fid: data = fid['data'] dts = fid['dts'] idx = np.logical_and(~np.isin(dts, dts0),dts>dt-timedelta(hours=3)) dts = np.hstack([dts[idx],dts0]) data = np.vstack([data[idx], data0]) idx = np.argsort(dts) dts = dts[idx] data = data[idx] else: dts = dts0 data = data0 np.savez(attfileout, data=data, dts=dts, vars=vars_list) def extract_orb(working_leo): orb_dir = working_dir / f'L1_extract_{working_leo}/rcvOrb' orbfiles = list(orb_dir.glob('*_sp3')) leo_sats = np.unique([sub(r'.*_\d{4}\.\d{3}\.(....)\.\d\d_.*', r'\1', x.name) for x in orbfiles]) # orbfiles = glob(f'/data3/xshao/planetiq/SWx/L1A/*/gnomes/level1a/rcvOrb/rcvOrb_*') # leo_sats = np.unique([sub(r'.*_\d{4}\.\d{3}\.(....)\.\d\d_.*', r'\1', x) for x in orbfiles]) # leo_dts = np.unique([datetime.strptime(sub(r'.*_(\d{4}\.\d{3})\..*',r'\1',x),'%Y.%j') for x in orbfiles]) # leo_sats = np.array(['GN04','GN05']) # leo_dts = np.array([datetime(2026, 1, 1) + timedelta(days=x) for x in range(47,102)]) for ff in leo_sats: log_message(f"-- Processing {ff}...", working_leo) orbfiles = list(orb_dir.glob(f'*{ff}*_sp3')) orbfiles.sort() dt = datetime.strptime(sub(r'.*_(\d{4}\.\d{3})\..*', r'\1', orbfiles[0].name), '%Y.%j') orbfileout = working_dir / f'p01_leoOrb/leoOrb_{ff}.{dt:%Y.%j}.npz' dts0 = [] data0 = [] tmp0 = np.full((8,), np.nan) for orbfile in orbfiles: for line in open(orbfile, 'r'): if line.startswith("*"): if dts0: data0.append(tmp0.copy()) tmp0[:] = np.nan dts0.append(to_datetime(line[1:].split())) elif line.startswith("P"): tmp0[:4] = np.array([line[x:x + 14] for x in range(4, 60, 14)]).astype('float') elif line.startswith("V"): tmp0[4:] = np.array([line[x:x + 14] for x in range(4, 60, 14)]).astype('float') data0.append(tmp0.copy()) dts0, idx = np.unique(np.array(dts0), return_index=True) data0 = np.array(data0)[idx] if orbfileout.exists(): with np.load(orbfileout, allow_pickle=True) as fid: data = fid['data'] dts = fid['dts'] idx = ~np.isin(dts, dts0) dts = np.hstack([dts[idx], dts0]) data = np.vstack([data[idx], data0]) idx = np.argsort(dts) dts = dts[idx] data = data[idx] else: dts = dts0 data = data0 np.savez(orbfileout, data=data, dts=dts) def gather_crx(working_leo): crx_dir = working_dir / f'L1_extract_{working_leo}/podRx3' from tools import get_sat_freq import pickle T0 = datetime(1980,1,6) header = ['GPS_seconds', 'time_offset', 'L1C', 'L2X', 'S1C', 'S2X', 'C1C', 'C2X'] crxfiles = list(crx_dir.glob('*_crx')) crxfiles.sort() # dtstr = sub(r'.*_(\d{4}\.\d{3})\.\d{3}\.\d\d\.\d\d_.*', r'\1', crxfiles[-1]) # dt = datetime.strptime(dtstr, '%Y.%j') # leo = lsat_names[int(sub(r'.*_\d{4}\.\d{3}\.(\d{3})\.\d\d\.\d\d_.*', r'\1', crxfiles[-1]))] # ant = int(sub(r'.*_\d{4}\.\d{3}\.\d{3}\.\d\d\.(\d\d)_.*', r'\1', crxfiles[-1])) dts = np.unique([datetime.strptime(sub(r'.*_(\d{4}\.\d{3})\.\d{3}\.\d\d\.\d\d_.*', r'\1', crxfile.name),'%Y.%j') for crxfile in crxfiles]) leos = np.unique([lsat_names[int(sub(r'.*_\d{4}\.\d{3}\.(\d{3})\.\d\d\.\d\d_.*', r'\1', crxfile.name))] for crxfile in crxfiles]) ants = np.unique([int(sub(r'.*_\d{4}\.\d{3}\.\d{3}\.\d\d\.(\d\d)_.*', r'\1', crxfile.name)) for crxfile in crxfiles]) # dtstrs = np.array([f'2026.{x:03d}' for x in range(48,103)]) # dts = np.array([datetime.strptime(dtstr0, '%Y.%j') for dtstr0 in dtstrs]) # leos = np.array(['GN04','GN05']) # ants = np.arange(2) for dt in dts: (working_dir / f'p01_podPair/{dt:%Y-%m-%d}').mkdir(exist_ok=True) for leo in leos: for ant in ants: crxfileout = working_dir / f'p01_podRx3/podRx3_{leo}.{dt:%Y.%j}.{ant:02d}.pkl' crxfileprev = working_dir / f'p01_podRx3/podRx3_{leo}.{dt-timedelta(days=1):%Y.%j}.{ant:02d}.pkl' if crxfileout.exists(): log_message(f'-- Combining {crxfileout.name}', working_leo) with open(crxfileout, 'rb') as fid: pkldata = pickle.load(fid) dsat = pkldata['dsat'] dtime = pkldata['dtime'] dataall = pkldata['dataall'] elif crxfileprev.exists(): log_message(f'-- Retrieving {crxfileprev.name}', working_leo) with open(crxfileprev, 'rb') as fid: pkldata = pickle.load(fid) dsat = pkldata['dsat'] dtime = pkldata['dtime'] dataall = pkldata['dataall'] for dnum in range(len(dsat)): tsidx = np.where(np.diff(np.hstack([[T0], dtime[dnum]])) > timedelta(seconds=3))[0].astype('int') teidx = np.hstack([tsidx[1:], tsidx[:1]]) - 1 tidx = tsidx[dtime[dnum][teidx] > dt - timedelta(hours=1)] if tidx.size > 0: dtime[dnum] = dtime[dnum][tidx[0]:] dataall[dnum] = dataall[dnum][tidx[0]:, :] else: # assign empty array with .size == 0 to be removed afterward dtime[dnum] = dtime[dnum][:0] dataall[dnum] = dataall[dnum][:0] didx = np.where([tt.size > 0 for tt in dtime])[0] dsat = [dsat[x] for x in didx] dtime = [np.array(dtime[x]) for x in didx] dataall = [np.array(dataall[x]) for x in didx] else: log_message(f'-- Creating {crxfileout.name}', working_leo) dsat = [] dtime = [] dataall = [] crxfiles = list(crx_dir.glob(f'podRx3_{dt:%Y.%j}.{lsat_names[leo]:03d}.*.{ant:02d}_crx')) crxfiles.sort() for crxfile in crxfiles: is_head = True for line in open(crxfile, 'r'): if is_head: if 'END OF HEADER' in line: is_head = False continue if line.startswith(">"): dt_now = to_datetime(line[1:].split()) elif line[0] in ('G','R','E','C'): gtmp = [line[ii:ii + 14] for ii in range(3, len(line) - 1, 16)] gtmp = [float(ss) if ss.strip() else np.nan for ss in gtmp] gtmp += [np.nan for _ in range(6 - len(gtmp))] gsat = line[:3] if gsat not in dsat: dsat.append(gsat) dtime.append(np.array([dt_now])) dataall.append(np.array([gtmp])) else: gnum = dsat.index(gsat) if dt_now > dtime[gnum][-1]: dtime[gnum] = np.hstack([dtime[gnum], [dt_now]]) dataall[gnum] = np.vstack([dataall[gnum],[gtmp]]) elif dt_now in dtime[gnum]: tidx = np.searchsorted(dtime[gnum], dt_now) dataall[gnum][tidx] = [pp if np.isfinite(pp) else qq for pp, qq in zip(dataall[gnum][tidx], gtmp)] else: tidx = np.searchsorted(dtime[gnum], dt_now) dtime[gnum] = np.hstack([dtime[gnum][:tidx], [dt_now], dtime[gnum][tidx:]]) dataall[gnum] = np.vstack([dataall[gnum][:tidx], [gtmp], dataall[gnum][tidx:]]) didx = np.argsort(dsat) dsat = [dsat[x] for x in didx] dtime = [np.array(dtime[x]) for x in didx] dataall = [np.array(dataall[x]) for x in didx] pkldata = {'dsat': dsat, 'dtime': dtime, 'dataall': dataall} with open(crxfileout, 'wb') as fid: pickle.dump(pkldata, fid) c0 = 299792458 for dnum in range(len(dsat)): f1,f2,_ = get_sat_freq(dsat[dnum]) idx = np.all(np.isfinite(dataall[dnum]), axis=1) dtime[dnum] = dtime[dnum][idx] dataall[dnum] = dataall[dnum][idx] idx = [0] + list(np.where(np.diff(dtime[dnum]) > timedelta(seconds=3))[0] + 1) + [len(dtime[dnum])] for ss in range(len(idx) - 1): if idx[ss + 1] - idx[ss] < 60: continue gps_sec = (dtime[dnum][idx[ss]] - T0) / timedelta(seconds=1) tdiff = (dtime[dnum][idx[ss]:idx[ss + 1]] - dtime[dnum][idx[ss]]) / timedelta(seconds=1) data = np.column_stack(( gps_sec + np.zeros_like(tdiff), tdiff, dataall[dnum][idx[ss]:idx[ss+1], [1, 4, 2, 5, 0, 3]] * np.array([c0 / f1, c0 / f2, 1, 1, 1, 1]) )).astype('float64') tStart, tEnd = [T0 + timedelta(seconds=x[0] + x[1]) for x in data[[0, -1], :2]] occFileName = working_dir / f'p01_podPair/{dt:%Y-%m-%d}/occTab_{leo}.{dt:%Y.%j}.{tStart:%Y%m%d%H%M%S}.{tEnd:%Y%m%d%H%M%S}.A{ant:02d}.{dsat[dnum]}.npz' if occFileName.exists(): continue for oldfile in (working_dir / f'p01_podPair/{dt:%Y-%m-%d}').glob(f'occTab_{leo}.{dt:%Y.%j}.{tStart:%Y%m%d%H%M%S}.*.A{ant:02d}.{dsat[dnum]}.npz'): token = sub(r'.*\.(\d{14})\.(\d{14})\.(A\d\d)\.([A-Z]\d\d)\.npz', r'\1_\2_\4_\3', oldfile.name) for partfile in (working_dir / f'p05_podTec/{dt:%Y-%m-%d}').glob(f'*{leo}_{token}*'): partfile.unlink() oldfile.unlink() np.savez(occFileName, data=data, header=header) def gather_data(working_leo): log_message('- Start gathering leoOrb profiles.', working_leo) extract_orb(working_leo) log_message('- Start gathering leoAtt profiles.', working_leo) extract_att(working_leo) log_message('- Start gathering REF profiles.', working_leo) gather_crx(working_leo)