from datetime import datetime, timedelta from glob import glob import gzip import os from subprocess import run from re import sub from gnss.tools import dt2float, lagrange #Jun modified 20260504 import numpy as np from gnss.tools import read_GNSS, gps_week #Jun modified 20260504 from configuration import gnss_dir,gnss_download_url1,gnss_download_url2 #Jun modified 20260504 # from matplotlib import pyplot as plt # from scipy.interpolate import CubicSpline # from frame_conversions import ro_ecef2eci, ro_eci2ecef # from astropy.time import Time # from astropy.coordinates import ITRS, GCRS, CartesianRepresentation # import astropy.units as u # Jun modified 20260504 begin -------------------------------- #os.chdir('/data3/tcliu/gps_pod/python/gnss') # original os.chdir(gnss_dir) os.makedirs(os.path.join(gnss_dir,'gnss5min'),exist_ok=True) os.makedirs(os.path.join(gnss_dir,'gnssLrg'),exist_ok=True) os.makedirs('log',exist_ok=True) # Jun modified 20260504 end ------------------------------- logfile = f"log/interp_gnss_{datetime.now():%Y%m%d}.txt" def log_message(msg): with open(logfile, 'a') as fid: fid.write(f"{datetime.now():%H:%M:%S} {msg}\n") log_message('- Downloading...') tar_dirs = glob('gnssLrg/*') # Jun modified 20260504: just in case we start from scratch and tar_dirs=[]------------ if len(tar_dirs)>0: tar_dirs.sort() tar_sdmin = datetime.strptime(sub(r'.*/',r'',tar_dirs[-1]),'%Y.%j')-timedelta(days=1) date_len = np.ceil((datetime.now() - tar_sdmin) / timedelta(days=1)).astype(int) + 1 date_rng = [tar_sdmin + timedelta(days=x) for x in range(date_len)] else: dt_today = datetime.today() date_rng = [dt_today + timedelta(days=x) for x in range(-5,2)] # Jun modified 20260504----------------------------------------- # dt_proc = datetime(2025, 11, 6) # dt_file = [dt_proc + timedelta(hours=x) for x in range(-120, 300, 6)] # dt_range = [dt_proc + timedelta(hours=x) for x in [-1,25]] # Jun modified 20260504 ---------------------------------------- #urls = np.hstack([[f'ftps://gdc.cddis.eosdis.nasa.gov/gnss/products/latest/ultra/GRG0OPSULT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz', # f'ftps://gdc.cddis.eosdis.nasa.gov/gnss/products/latest/ultra/COD0OPSULT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz', # f'ftps://gdc.cddis.eosdis.nasa.gov/gnss/products/latest/ultra/GFZ0OPSULT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz', # f'ftps://gdc.cddis.eosdis.nasa.gov/gnss/products/{gps_week(x):d}/WUM0MGXNRT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz'] for x in date_rng]) urls = np.hstack([[f'{gnss_download_url1}latest/ultra/GRG0OPSULT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz', f'{gnss_download_url1}latest/ultra/COD0OPSULT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz', f'{gnss_download_url1}latest/ultra/GFZ0OPSULT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz', f'{gnss_download_url1}{gps_week(x):d}/WUM0MGXNRT_{x:%Y%j}*00_02D_05M_ORB.SP3.gz'] for x in date_rng]) #Jun modified 20260504 ----------------------------------------------- # urls = ['ftps://gdc.cddis.eosdis.nasa.gov/gnss/products/latest/final/GRG0OPSFIN_20253*00_01D_30S_CLK.CLK.gz', # 'ftps://gdc.cddis.eosdis.nasa.gov/gnss/products/latest/final/GRG0OPSFIN_20253*00_01D_05M_ORB.SP3.gz', # 'ftps://gdc.cddis.eosdis.nasa.gov/gnss/products/latest/final/COD0OPSFIN_20253*00_01D_05M_ORB.SP3.gz'] for ftps_url in urls: command = [ "wget", "-N", "--ftp-user=anonymous", # Specify FTP username "--ftp-password=passwd", # Specify FTP password ftps_url ] run(command, cwd='./gnss5min/', capture_output=True) log_message('- Download is completed.') # Jun modified 20260504 -------------------------------- #urls = [f'https://gpsmet.umd.edu/GNSSOrb_Ext/GNSSOrb_{x-timedelta(days=2):%Y%j}0_withCLK.SP3' for x in date_rng] #original urls = [f'{gnss_download_url2}GNSSOrb_{x-timedelta(days=2):%Y%j}0_withCLK.SP3' for x in date_rng] # Jun modified 20260504 ---------------------------- for ftps_url in urls: command = ['wget',ftps_url] run(command, cwd='./gnss5min/', capture_output=True) lfile = sub(r'.*/GNSSOrb_(\d{7})0_.*', r'GBM0MGXRAP_\g<1>0000_04D_05M_ORB.SP3', ftps_url) command = ['mv', sub(r'.*/', r'', ftps_url), lfile] run(command, cwd='./gnss5min/', capture_output=True) command = ['gzip', lfile] run(command, cwd='./gnss5min/', capture_output=True) log_message('- Download from Bin folder on gpsmet is completed.') # # for dt_roll in dt_file: # # https_url = f'https://isdc-data.gfz.de/gnss/products/ultra/w{(dt_roll-datetime(1980,1,6))//timedelta(weeks=1)}/GFZ0OPSULT_{dt_roll:%Y%j%H}00_02D_05M_ORB.SP3.gz' # # run(["wget","-N",https_url], cwd='./gnss5min/', capture_output=True) # https_url = f'https://isdc-data.gfz.de/gnss/products/final/w{(dt_roll - datetime(1980, 1, 6)) // timedelta(weeks=1)}/GFZ0OPSFIN_{dt_roll:%Y%j%H}00_01D_15M_ORB.SP3.gz' # run(["wget", "-N", https_url], cwd='./gnss5min/', capture_output=True) # https_url = f'https://isdc-data.gfz.de/gnss/products/final/w{(dt_roll - datetime(1980, 1, 6)) // timedelta(weeks=1)}/GFZ0OPSFIN_{dt_roll:%Y%j%H}00_01D_30S_CLK.CLK.gz' # run(["wget", "-N", https_url], cwd='./gnss5min/', capture_output=True) gdir = './gnss5min/' idir = './gnssLrg/' gfiles = glob(gdir + '*SP3.gz') gfiles.sort() for gfile in gfiles: dt = datetime.strptime(sub(r'.*_(\d{7})\d{4}_.*',r'\1',gfile),'%Y%j') tdir = f'{gdir}{dt:%Y.%j}' os.makedirs(tdir, exist_ok=True) ldir = f'{idir}{dt:%Y.%j}' os.makedirs(ldir, exist_ok=True) tfile = sub(r'.*/(.*)\.gz',f'{tdir}/\\1.txt',gfile) ftoken = sub(r'.*/(.*)0[A-Z]{6}(_\d{9}).*', r'\1\2', tfile) hr = int(ftoken[-2:]) if hr%6 != 0: os.remove(gfile) continue ofile = glob(f'{ldir}/{ftoken}_*.npz') mdt = datetime.fromtimestamp(os.path.getmtime(gfile)) if len(ofile) > 0: odt = datetime.strptime(sub(r'.*_c(\d+)\.npz',r'\1',ofile[-1]), '%Y%j%H%M') if mdt - odt > timedelta(minutes=2): os.remove(ofile[-1]) log_message(f'- Updating {ftoken} ({odt:%m-%d %H:%M} >> {mdt:%H:%M})...') else: os.remove(gfile) log_message(f'- Skipping {ftoken} ({odt:%m-%d %H:%M})...') continue else: log_message(f'- Creating {ftoken} ({mdt:%m-%d %H:%M})...') with gzip.open(gfile,'rb') as fin, open(tfile,'wb') as fout: fout.write(fin.read()) [Gdata0, Gdt0, Gsats] = read_GNSS(tfile) sidx = np.all(np.isfinite(Gdata0[:,:,:4]),axis=(0,2)) Gsats = Gsats[sidx] Gdata0 = Gdata0[:,sidx,:] tidx = np.argsort(Gdt0) sidx = np.argsort(Gsats) Gdt0 = Gdt0[tidx] Gsats = Gsats[sidx] Gdata0 = Gdata0[np.ix_(tidx, sidx)] g2file = sub(r'ORB\.SP3',r'CLK.CLK',gfile) if os.path.isfile(g2file): t2file = sub(r'.*/(.*)\.gz',f'{tdir}/\\1.txt',g2file) with gzip.open(g2file, 'rb') as fin, open(t2file, 'wb') as fout: fout.write(fin.read()) for line in open(t2file,'r'): if line.startswith('AS'): tt = datetime.strptime(line[8:24],'%Y %m %d %H %M') ss = line[3:6] tidx = np.searchsorted(Gdt0, tt) sidx = np.searchsorted(Gsats, ss) if Gdt0[tidx] == tt and Gsats[sidx] == ss: Gdata0[tidx,sidx,3] = float(line[40:59])*1E6 Gdt = np.arange(Gdt0[0],Gdt0[-1]+timedelta(seconds=30),timedelta(seconds=30)).astype(datetime) p, v = lagrange(dt2float(Gdt0, dt) / 86400, Gdata0[:, :, :4].reshape((Gdata0.shape[0], -1)), dt2float(Gdt, dt) / 86400, 11) Gdata = np.concatenate((p.reshape(p.shape[0], -1, 4), (v * 10000 / 86400).reshape(v.shape[0], -1, 4)), axis=2) ofile = f'{ldir}/{ftoken}_s{Gdt[0]:%Y%j%H%M}_e{Gdt[-1]:%Y%j%H%M}_c{mdt:%Y%j%H%M}.npz' np.savez(ofile, Gdata=Gdata, Gdt=Gdt, Gsats=Gsats) os.remove(tfile) os.replace(gfile, tfile[:-3]+'gz') if os.path.isfile(g2file): os.remove(t2file) os.replace(g2file, t2file[:-3] + 'gz') log_message('- Progress is completed.') # gdir = './gnss5min/' # idir = './gnssFin/' # gfiles = glob(gdir + '*SP3.gz') # for gfile in gfiles: # dt = datetime.strptime(sub(r'.*_(\d{7})\d{4}_.*',r'\1',gfile),'%Y%j') # tdir = f'{gdir}{dt:%Y.%j}' # os.makedirs(tdir, exist_ok=True) # ldir = f'{idir}{dt:%Y.%j}' # os.makedirs(ldir, exist_ok=True) # tfile = sub(r'.*/(.*)\.gz',f'{tdir}/\\1.txt',gfile) # with gzip.open(gfile,'rb') as fin, open(tfile,'wb') as fout: # fout.write(fin.read()) # [Gdata, Gdt, Gsats] = read_GNSS(tfile) # ofile = sub(r'.*/(.*)\.gz',f'{ldir}/\\1.npz',gfile) # np.savez(ofile, Gdata=Gdata, Gdt=Gdt, Gsats=Gsats) # os.remove(tfile) # os.replace(gfile, tfile[:-3]+'gz') # for ss in ['COD','GRG','GFZ']: # Gdatam = np.empty((0,0,8)) # Gdtm = np.empty((0,)) # Gsats = np.empty((0,)) # Gdt = np.array([dt_proc + timedelta(seconds=x) for x in range(0,86430,30)]) # for dt_roll in dt_file: # # dt_now = None # gfile = f'./gnss5min/{dt_roll:%Y.%j}/{ss}0OPSULT_{dt_roll:%Y%j%H}00_02D_05M_ORB.SP3.txt' # if not os.path.exists(gfile): # continue # [Gdata0, Gdt0, Gsats0] = read_GNSS(gfile) # tidx = np.logical_and(Gdt0 >= dt_range[0], Gdt0 <= dt_range[1]) # Gdata0 = Gdata0[tidx] # Gdt0 = Gdt0[tidx] # tidx = np.argsort(Gdt0) # sidx = np.argsort(Gsats0) # Gdt0 = Gdt0[tidx] # Gsats0 = Gsats0[sidx] # Gdata0 = Gdata0[np.ix_(tidx,sidx)] # sidx = ~np.isin(Gsats0, Gsats) # Gsats = np.hstack((Gsats, Gsats0[sidx])) # Gdatam = np.hstack((Gdatam, np.full((Gdatam.shape[0],np.sum(sidx),Gdatam.shape[2]),np.nan))) # tidx = ~np.isin(Gdt0, Gdtm) # Gdtm = np.hstack((Gdtm, Gdt0[tidx])) # Gdatam = np.vstack((Gdatam, np.full((np.sum(tidx),Gdatam.shape[1],Gdatam.shape[2]),np.nan))) # tidx = np.argsort(Gdtm) # sidx = np.argsort(Gsats) # Gdtm = Gdtm[tidx] # Gsats = Gsats[sidx] # Gdatam = Gdatam[np.ix_(tidx,sidx)] # Gdatam[np.ix_(np.isin(Gdtm, Gdt0),np.isin(Gsats, Gsats0))] = Gdata0 # # tidx = Gdt