import sys,os import netCDF4,glob,collections,pickle import numpy as np import matplotlib #matplotlib.use('Agg') #matplotlib.use('TkAgg') import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cf from matplotlib.colors import LogNorm from matplotlib.ticker import LogLocator, NullFormatter from matplotlib.ticker import FixedLocator class plot_BA_Tdry_class: def __init__(self,mission,leoid, ymd, input_dir, output_dir): self.mission,self.leoid,self.ymd,self.input_dir,self.output_dir=mission,leoid, ymd, input_dir, output_dir self.H_ba=[5,10,20,30] self.H_tdry=[15,20,25,30] self.out_ba={h:[] for h in self.H_ba} self.out_tdry={h:[] for h in self.H_tdry} def ingest(self,read=True): if read: if self.leoid=='all': fs=sorted(glob.glob(os.path.join(self.input_dir,self.ymd,'*.nc'))) else: fs=sorted(glob.glob(os.path.join(self.input_dir,self.ymd,'RONA_L1C_'+self.leoid+'*.nc'))) recs=[] for f in fs: d=netCDF4.Dataset(f) d.set_auto_mask(False) bad=str(d.getncattr('bad')).strip() if bad!='0': d.close() continue recs.append(dict(f=f,lat=float(d['lat'][0]), #reference lon=float(d['lon'][0]), roc=float(d['roc'][0]), und=float(d['undulation'][0]), ip=d['impact_opt'][0].copy(), #impact parameter (m) >0 dim_lev1b ba=d['bangle_opt'][0].copy(), #baneding angle (radians) >0 dim_lev1b lat_tp=d['lat_tp'][0].copy(), # lat at TP (degree) -90~90 dim_lev1b lon_tp=d['lon_tp'][0].copy(), # lon at TP (degree) -180~180 dim_lev1b z=d['alt_refrac'][0].copy(), # geomeric height above geoid (m) -1000~150000 dim_lev2a tdry=d['dry_temp'][0].copy())) # Tdry (K) 150~500 dim_lev2a d.close() pickle.dump(recs,open(os.path.join(self.output_dir,f'recs_{self.mission}_{self.leoid}_{self.ymd}.pkl'),'wb')) else: recs=pickle.load(open(os.path.join(self.output_dir,f'recs_{self.mission}_{self.leoid}_{self.ymd}.pkl'),'rb')) LA,LO,BA=[],[],[] # for global all points for BA for irec,r in enumerate(recs): # BA on impact height ========================= v=(r['ba']>=0)&(r['ba']<=0.1)&(r['ip']>=6.2e6)&(r['ip']<=6.6e6)&(r['lat_tp']>=-90)&(r['lat_tp']<=90)&(r['lon_tp']>=-180)&(r['lon_tp']<=180) if v.sum()<2: continue ih=(r['ip']-r['roc']-r['und'])/1e3 # impact height km ih,ba,la,lo=ih[v],r['ba'][v],r['lat_tp'][v],r['lon_tp'][v] LA.append(la);LO.append(lo);BA.append(ba) # save for global all points # increase, for interpolation o=np.argsort(ih) ih,ba,la,lo=ih[o],ba[o],la[o],lo[o] lou=np.degrees(np.unwrap(np.radians(lo))) # lon unpacked for h in self.H_ba: if ih[0]<=h<=ih[-1]: self.out_ba[h].append((np.interp(h,ih,la),(np.interp(h,ih,lou)+180)%360-180,np.exp(np.interp(h,ih,np.log(ba))))) # Tdry on alt_refac=========================== z=r['z'] #km t=r['tdry'] #K v=(z>=-1000)&(z<=150000)&(t>=150)&(t<=500) if v.sum()<2: continue z,t=z[v]/1e3,t[v] o=np.argsort(z) z,t=z[o],t[o] #print(irec,r['f'],z) for h in self.H_tdry: if z[0]<=h<=z[-1]: T=np.interp(h,z,t) y,x=np.interp(h,ih,la),(np.interp(h,ih,lou)+180)%360-180 self.out_tdry[h].append((y,x,T)) self.LA,self.LO,self.BA=map(np.concatenate,(LA,LO,BA)) self.nprf=len(recs) print('------') def plots(self): #self.plot_allalts(task='BA') self.plot_4levels(task='BA') #self.plot_4levels(task='Tdry') def base_old(self,ax): ax.set_global(); ax.add_feature(cf.LAND.with_scale('50m'),facecolor='#e8e8e8',zorder=0); ax.coastlines('50m',lw=0.5,color='0.3') #gl=ax.gridlines(draw_labels=True,lw=0.3,color='0.6',alpha=0.6) gl=ax.gridlines(draw_labels={'bottom': 'x','left': 'y'},lw=0.3,color='0.6',alpha=0.6) # Latitude every 30 degrees gl.ylocator = FixedLocator([-90, -60, -30, 0, 30, 60, 90]) gl.xlocator = FixedLocator([-180, -120, -60, 0, 60, 120, 180]) gl.top_labels = False gl.right_labels = False gl.left_labels = True #gl.geo_labels = False gl.xlabel_style=gl.ylabel_style={'size':8} def base(self, ax): ax.set_global() ax.add_feature( cf.LAND.with_scale('50m'), facecolor='#e8e8e8', zorder=0 ) ax.coastlines('50m', lw=0.5, color='0.3') gl = ax.gridlines( draw_labels=True, lw=0.3, color='0.6', alpha=0.6 ) # Grid lines gl.xlocator = FixedLocator([-180, -120, -60, 0, 60, 120, 180]) gl.ylocator = FixedLocator([-60, -30, 0, 30, 60]) # Longitude labels at bottom gl.top_labels = False gl.bottom_labels = True # Turn off Cartopy's latitude labels gl.left_labels = False gl.right_labels = False gl.geo_labels = False gl.xlabel_style = {'size': 10} # Manually put latitude labels on LEFT side for lat in [-60, -30, 0, 30, 60]: x, y = ax.projection.transform_point( -180, lat, ccrs.PlateCarree() ) if lat > 0: label = f'{lat}°N' elif lat < 0: label = f'{abs(lat)}°S' else: label = '0°' ax.annotate( label, xy=(x, y), xytext=(-5, 0), textcoords='offset points', ha='right', va='center', fontsize=10 ) def plot_allalts(self,task='BA'): if task=='BA': x=self.LO y=self.LA z=self.BA*1e3 #print(z.min(),z.max()) title='Bending Angle (mrad)' crange=[1e-8,1e1] # Fig 1: all tangent points o=np.argsort(z) fig=plt.figure(figsize=(14,7.5)) plt.rcParams.update({'font.size': 12}) ax=plt.axes(projection=ccrs.Robinson()) self.base(ax) sc=ax.scatter(x[o],y[o],c=z[o],s=0.3,cmap='turbo',norm=LogNorm(crange[0],crange[1]),transform=ccrs.PlateCarree(),rasterized=True,lw=0) cb=plt.colorbar(sc,ax=ax,orientation='horizontal',pad=0.06,fraction=0.04,aspect=50,extend='both') cb.set_label('Bending angle, bangle_opt (mrad, log scale)') # Add minor ticks at 2,3,...,9 × 10^n # Minor log ticks: 2, 3, ..., 9 × 10^n cb.minorlocator = LogLocator(base=10,subs=np.arange(2, 10)) cb.update_ticks() cb.ax.xaxis.set_major_locator(LogLocator(base=10, subs=(1.0,), numticks=100)) cb.ax.xaxis.set_minor_locator(LogLocator(base=10, subs=np.arange(2, 10), numticks=100)) cb.ax.xaxis.set_minor_formatter(NullFormatter()) cb.ax.tick_params(axis='x',which='minor',bottom=True,top=False,length=3,width=1) cb.ax.tick_params(axis='x',which='major',bottom=True,top=False,length=6,width=1) ax.set_title(f'{self.mission} {title} {self.ymd}\n {self.nprf} profiles, {z.size:,} points',fontsize=18,y=1.02) outfig=os.path.join(self.output_dir,f'{self.mission}_{task}_{self.ymd}_allpt.png') plt.savefig(outfig,dpi=170,bbox_inches='tight') plt.show() plt.close() def plot_4levels(self,task='Tdry'): if task=='BA': out=self.out_ba cb_label='Bending Angle (mrad)' title='Bending Angle (mrad)' hname='Impact height' H=self.H_ba scale=1e3 cranges=[[12,19],[7,8],[1,2],[0.1,0.4]] #cranges=[[10,20],[7,8],[1,2],[0.3,0.4]] elif task=='Tdry': out=self.out_tdry cb_label='Dry temperature (K)' title='Dry temperature (K)' hname='MSL height' H=self.H_tdry scale=1 cranges=[[195,225],[195,225],[205,230],[215,240]] fig,axs=plt.subplots(2,2,figsize=(16,9.5),subplot_kw=dict(projection=ccrs.Robinson())) for ax,h,crange in zip(axs.flat,H,cranges): self.base(ax) a=np.array(out[h]) vmin,vmax=np.percentile(a[:,2],[1,99]) vmin=crange[0] vmax=crange[1] print(np.min(a[:,2]*scale),np.max(a[:,2]*scale)) sc=ax.scatter(a[:,1],a[:,0],c=a[:,2]*scale,s=7,cmap='turbo',vmin=vmin,vmax=vmax,transform=ccrs.PlateCarree(),edgecolor='k',lw=0.1) cb=plt.colorbar(sc,ax=ax,orientation='horizontal',pad=0.07,fraction=0.045,aspect=40,extend='both'); cb.set_label(cb_label,fontsize=10) ax.set_title(f'{hname} {h} km (n = {len(a)})',fontsize=14) fig.suptitle(f'{self.mission} {title} {self.ymd}',fontsize=16,y=0.96) plt.subplots_adjust(wspace=0.1, # horizontal space between panels hspace=0.3) # vertical space between panels plt.savefig(os.path.join(self.output_dir,f'{self.mission}_{task}_{self.ymd}_4levels.png'),dpi=170,bbox_inches='tight') #[mission,leoid, ymd, input_dir, output_dir]=sys.argv[1:6] mission,leoid, ymd='COSMIC2','all','2026-10-03' #mission,leoid, ymd='PlanetiQ','YM08','2026-09-29' input_dir=os.path.join(r'C:\Users\jzhou128\Desktop\workspace\GPS-RO\Topics\RO_pipeline\plots_nccf',mission.lower()+'rt','atmPrf_nccf') output_dir=input_dir print(mission,leoid, ymd, input_dir, output_dir) plot_BA_Tdry_obj=plot_BA_Tdry_class(mission,leoid, ymd, input_dir, output_dir) plot_BA_Tdry_obj.ingest(read=False) plot_BA_Tdry_obj.plots()