import os from pathlib import Path # import subprocess # from glob import glob from re import sub import numpy as np from datetime import datetime, timedelta from cosmic2.tools import to_datetime, get_sat_freq from struct import unpack, calcsize working_dir = Path(__file__).absolute().parent # lsat_names = dict(C2E1=1, C2E2=2, C2E3=3, C2E4=4, C2E5=5, C2E6=6) # for keys in list(lsat_names.keys()): # lsat_names[lsat_names[keys]] = keys def log_message(msg, working_leo): with open(working_dir / 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 gather_data(working_leo, tarfiles, filetypes, leo): log_message(f'- Start gathering RO profiles of {leo}.', working_leo) obs_type = 'CLDSXXXX' T0 = datetime(1980, 1, 6) gnsfiles = tarfiles[filetypes==3] log_message(f"-- {len(gnsfiles)} files to process", working_leo) tStart_min = datetime.now() + timedelta(days=1) tEnd_max = T0 for gnsfile in gnsfiles: dt = datetime.strptime(sub(r'.*_(\d{4}\.\d{3})\.\d\d\.\d\d\.\d{3}\.\d\d\..*',r'\1',gnsfile.name),'%Y.%j') # ant = int(sub(r'.*_\d{4}\.\d{3}\.\d\d\.\d\d\.\d{3}\.(\d\d)\..*',r'\1',gnsfile.name)) nRO = 0 with open(gnsfile, 'rb') as fid: fid.seek(-9,2) ver = unpack('>B', fid.read(1))[0] sats = unpack('>8s', fid.read(8))[0].strip().decode() if ver == 2: sat_preset = 36 elif ver == 3: sat_preset = 64 hr_offset = dict() end_byte = fid.seek(-(sat_preset * 8 + 4) * len(sats) - 9, 2) start_byte = {end_byte: 'EOF'} for sat in sats: hr_offset[sat] = unpack('>i', fid.read(4))[0] / 1E9 token = unpack(f'>{sat_preset}q', fid.read(sat_preset * 8)) for ii in range(len(token)): if token[ii] != -1: start_byte[token[ii]] = f'{sat}{ii+1:02d}' # no need, there are one-or-more profiles per gnss fid.seek(0) new_record = True while fid.tell() < end_byte: gsec, status, ant, gtype, svn, prn, rate, nlr, nhr = unpack('>IBBsHBHBB', fid.read(14)) fid.seek(nlr * 2, 1) # low rate header (skip) if new_record: header = ["GPS_seconds", "time_offset"] hr_dtype = '>' for _ in range(nhr): bits = unpack('>H',fid.read(2))[0] header.append(f"{obs_type[bits>>13]}{(bits>>10)%8+1}{64+(bits>>5)%32:c}"+ (f" ({64+bits%32:c})" if bits%32 else "")) hr_dtype += fid.read(1).decode() hr_dtype = sub('S', 'H', sub('s', 'h', hr_dtype)) # 2-byte int hr_dlen = calcsize(hr_dtype) data = [] gsec_0 = gsec gsat = f'{gtype.decode()}{prn:02d}' else: fid.seek(nhr * 3, 1) gsec_diff = gsec - gsec_0 + hr_offset[gtype.decode()] fid.seek(nlr * 8, 1) # low rate data (skip) for _ in range(rate): data.append((gsec_0, int.from_bytes(fid.read(3), 'big') / 1E7 + gsec_diff) + unpack(hr_dtype, fid.read(hr_dlen))) gsec_n, status_n, _, gtype_n, _, prn_n, _, _, _ = unpack('>IBBsHBHBB', fid.read(14)) fid.seek(-14, 1) # if fid.tell() in start_byte or gsec_n-gsec<1 or gsec_n-gsec>5 or gtype_n!=gtype or prn_n!=prn or status_n==0: if fid.tell() == end_byte or not np.logical_and(gsec_n-gsec==1,prn_n==prn): # TC: updated 2026/07/27 (add 1 line and indent 9 following lines) if 'L2P' not in header: data = np.array(data) tStart, tEnd = [T0 + timedelta(seconds=gsec_0 + data[x,1]) for x in [0, -1]] if tStart < tStart_min: tStart_min = tStart if tEnd > tEnd_max: tEnd_max = tEnd occFileName = working_dir / f"t01_RO_{working_leo}/occTab_{leo}.{dt:%Y.%j}.{tStart:%Y%m%d%H%M%S}.{tEnd:%Y%m%d%H%M%S}.A{ant:02d}.{gsat}.npz" np.savez(occFileName, data=data, header=header, source=gnsfile) nRO += 1 new_record = True else: new_record = False log_message(f"-- Processed {gnsfile.name}: {nRO} profiles", working_leo) log_message('- Start gathering REF profiles.', working_leo) import pickle c0 = 299792458 header = ['GPS_seconds', 'time_offset', 'L1C', 'L2X', 'S1C', 'S2X', 'C1C', 'C2X'] crxfiles = tarfiles[filetypes==2] # 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.array([datetime.strptime(crxfile.parent.as_posix()[-8:], '%Y/%j') for crxfile in crxfiles]) uniq_dts = np.unique(dts) # 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.array([int(sub(r'.*_\d{4}\.\d{3}\.\d\d\.\d\d\.\d{3}\.(\d\d)\..*', r'\1', crxfile.name)) for crxfile in crxfiles]) # ants = np.array([int(sub(r'.*_\d{4}\.\d{3}\.\d{3}\.\d\d\.(\d\d)_.*', r'\1', crxfile.name)) for crxfile in crxfiles]) uniq_ants = np.unique(ants) for dt in uniq_dts: # for leo in leos: for ant in uniq_ants: crxfileout = working_dir / f's01_podRx3/podRx3_{leo}.{dt:%Y.%j}.{ant:02d}.pkl' crxfileprev = working_dir / f's01_podRx3/podRx3_{leo}.{dt-timedelta(days=1):%Y.%j}.{ant:02d}.pkl' if os.path.exists(crxfileout): log_message(f'-- Combining {crxfileout}', working_leo) with open(crxfileout, 'rb') as fid: pkldata = pickle.load(fid) dsat = pkldata['dsat'] dtime = pkldata['dtime'] dataall = pkldata['dataall'] elif os.path.exists(crxfileprev): log_message(f'-- Retrieving {crxfileprev}', 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] > min(tStart_min, 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}', working_leo) dsat = [] dtime = [] dataall = [] crxfiled = crxfiles[np.logical_and(dts == dt, ants == ant)] crxfiled.sort() for crxfile in crxfiled: log_message(f'--- Processing {crxfile}', working_leo) is_head = True count = 0 sat_append = False for line in open(crxfile, 'r'): if is_head: if 'END OF HEADER' in line: is_head = False continue if sat_append: sats = [line[32+3*x:35+3*x] for x in range(min(12, sat_num))][::-1] + sats if sat_num>12: sat_num -= 12 else: sat_append = False elif count == 0: dt_now = to_datetime(('20'+line[1:]).split()) sat_num = int(line[29:32]) sats = [line[32+3*x:35+3*x] for x in range(min(12,sat_num))][::-1] count = sat_num*2 if sat_num>12: sat_num -= 12 sat_append = True elif count%2==0: count -= 1 gtmp = [line[ii:ii + 14] for ii in range(0, len(line) - 1, 16)] else: count -= 1 gtmp += [line[ii:ii + 14] for ii in range(0, len(line) - 1, 16)] gtmp = [float(ss) for ss in gtmp if ss.strip()] gsat = sats.pop() if gsat.startswith('G') or gsat.startswith('E'): 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} log_message(f'--- Saving {crxfileout}', working_leo) with open(crxfileout, 'wb') as fid: pickle.dump(pkldata, fid) rofiles = list((working_dir / f't01_RO_{working_leo}').glob('*.npz')) dts = np.array([datetime.strptime(sub(r'.*\.(\d{4}\.\d{3})\..*',r'\1',rofile.name), '%Y.%j') for rofile in rofiles]) uniq_dts = np.unique(dts) for dt in uniq_dts: crxfiles = list((working_dir / 's01_podRx3').glob(f'podRx3_{leo}.{dt:%Y.%j}*.pkl')) for crxfile in crxfiles: ant = int(sub(r'.*\.(\d\d)\.pkl',r'\1',crxfile.name)) with open(crxfile, 'rb') as fid: pkldata = pickle.load(fid) dsat = pkldata['dsat'] dtime = pkldata['dtime'] dataall = pkldata['dataall'] 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 or dtime[dnum][idx[ss + 1] - 1] < tStart_min or dtime[dnum][idx[ss]] > tEnd_max: 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], [0, 1, 4, 5, 2, 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't01_REF_{working_leo}/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' np.savez(occFileName, data=data, header=header) log_message('- Start gathering leoOrb profiles.', working_leo) orbfiles = tarfiles[filetypes==0] orbfiles.sort() dts_orb = np.array([datetime.strptime(sub(r'(\d{4}).(\d{3})',r'\1/\2',orbfile.parent.as_posix()[-8:]), '%Y/%j') for orbfile in orbfiles]) uniq_dts = np.unique(dts_orb) for dt in uniq_dts: log_message(f"-- Processing {dt:%Y-%m-%d}...", working_leo) orbfileout = working_dir / f's01_leoOrb/leoOrb_{leo}.{dt:%Y.%j}.npz' orbfileprev = working_dir / f's01_leoOrb/leoOrb_{leo}.{dt-timedelta(days=1):%Y.%j}.npz' orbfiled = orbfiles[dts_orb==dt] dts0 = [] data0 = [] tmp0 = np.full((8,),np.nan) for orbfile in orbfiled: isheader = True for line in open(orbfile,'r'): if line.startswith("*"): if not isheader: data0.append(tmp0.copy()) else: isheader = False 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 = np.array(dts0) data0 = np.array(data0) idx = np.abs(data0[:, 3] - 1E6) > 1E-3 dts0 = dts0[idx] data0 = data0[idx] # TC: updated 2027/07/27 (modified 1 line into 3 lines) # dts0, idx = np.unique(np.array(dts0), return_index=True) _, idx = np.unique(dts0[::-1], return_index=True) idx = len(dts0) - 1 - idx dts0 = dts0[idx] data0 = np.array(data0)[idx] if orbfileout.is_file(): 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] elif orbfileprev.is_file(): with np.load(orbfileprev, 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(orbfileout, data=data, dts=dts) log_message('- Start gathering leoAtt profiles.', working_leo) attfiles = tarfiles[filetypes==1] vars_list = ['att_x', 'att_y', 'att_z', 'att_w', 'sca_x', 'sca_y', 'sca_z', 'sca_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', 'tor', 'sad_x','sad_y'] dts_att = np.array([datetime.strptime(attfile.parent.as_posix()[-8:], '%Y/%j') for attfile in attfiles]) uniq_dts = np.unique(dts_att) for dt in uniq_dts: log_message(f"-- Processing {dt:%Y-%m-%d}...", working_leo) attfiled = attfiles[dts_att==dt] attfiled.sort() attfileout = working_dir / f'./s01_leoAtt/leoAtt_{leo}.{dt:%Y.%j}.npz' attfileprev = working_dir / f'./s01_leoAtt/leoAtt_{leo}.{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("sca"): tmp0[4:8] = np.fromstring(line[3:], sep=" ") elif line.startswith("ang"): tmp0[8:11] = np.fromstring(line[3:], sep=" ") elif line.startswith("pve"): tmp0[11:17] = np.fromstring(line[3:], sep=" ") elif line.startswith("pvi"): tmp0[17:23] = np.fromstring(line[3:], sep=" ") elif line.startswith("tor"): tmp0[23:24] = np.fromstring(line[3:], sep=" ") elif line.startswith("sad"): tmp0[24:] = np.fromstring(line[3:], sep=" ") data0.append(tmp0.copy()) dts0 = np.array(dts0) data0 = np.array(data0) if os.path.isfile(attfileout): 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 os.path.isfile(attfileprev): 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) log_message('- Completed gathering profiles.', working_leo)