#!/usr/bin/env python #Plot output from retrieval for whymsie flights. #Discontinuities in initial TB error, suspect pseudoindex interp is switched around. #Try plotting t field or sst across 25.75 lat #add clw import pdb import sys import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import numpy as np import matplotlib.ticker as ticker def myplot(lats,lons,var,title,bartitle,vmin,vmax,extent,dotsize=8,date=''): # 3. Create the Map Plot fig = plt.figure(figsize=(6, 8)) # Use the PlateCarree projection (standard lat/lon) ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) # Add geographic features for context ax.add_feature(cfeature.LAND, facecolor='lightgray', alpha=0.5) ax.add_feature(cfeature.OCEAN, facecolor='lightblue', alpha=0.5) ax.add_feature(cfeature.COASTLINE, linewidth=1) # Create a scatter plot of the data # We use c=statuses to color-code by the status variable #print('##',len(lons),len(var) ) scatter = ax.scatter(lons, lats, c=var, cmap='viridis', transform=ccrs.PlateCarree(), s=dotsize, edgecolor='none', zorder=5, vmin=vmin,vmax=vmax) fig.suptitle(title+'\n'+date,color='black') # Add a colorbar to interpret the status colors plt.colorbar(scatter, label=bartitle, shrink=0.5) # Set the map extent to frame the data properly with a little padding # Looking at the coordinates, this data is located off the coast of Southern California padding = 0.05 #ax.set_extent([ min(lons) - padding, max(lons) + padding, # min(lats) - padding, max(lats) + padding ], crs=ccrs.PlateCarree()) ax.set_extent(extent) # Add gridlines and coordinate labels gl = ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False, linestyle='--', color='gray', alpha=0.7) gl.top_labels = False gl.right_labels = False plt.title(title) plt.ion() plt.show() return(fig) def plot_side_by_side_maps(lat, lon, var1, var2, title1="Variable 1", title2="Variable 2", cbar_label="Data Values", cmap='viridis'): #From ChatGSFC """ Plots two variables side-by-side on PlateCarree maps with a corrected aspect ratio and a common colorbar on the right. Parameters: lat, lon : array-like coordinates var1, var2 : array-like data variables corresponding to lat/lon title1, title2 : str, titles for the respective subplots cbar_label : str, label for the common colorbar cmap : str, colormap name """ # 1. Calculate common min and max to ensure colors match perfectly across both maps vmin = min(np.nanmin(var1), np.nanmin(var2)) vmax = max(np.nanmax(var1), np.nanmax(var2)) # 2. Calculate the geographically accurate aspect ratio based on the center latitude # Aspect = 1 / cos(latitude). We use absolute value to handle southern hemisphere properly. center_lat = np.nanmean(lat) aspect_ratio = 1.0 / np.cos(np.radians(abs(center_lat))) # 3. Initialize Figure using Constrained Layout to protect the titles fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6), subplot_kw={'projection': ccrs.PlateCarree()}, layout='constrained') # 4. Plot Variable 1 # Note: lon is X, lat is Y. print(np.shape(lon)) print(np.shape(var1)) print(np.shape(var2)) sc1 = ax1.scatter(lon, lat, c=var1, cmap=cmap, vmin=vmin, vmax=vmax, transform=ccrs.PlateCarree(), s=15 ) #,alpha=0.9) ax1.coastlines() ax1.set_title(title1,y=1.02) ax1.set_aspect(aspect_ratio) # Apply the corrected aspect ratio # 5. Plot Variable 2 sc2 = ax2.scatter(lon, lat, c=var2, cmap=cmap, vmin=vmin, vmax=vmax, transform=ccrs.PlateCarree(), s=15) ax2.coastlines() ax2.set_title(title2,y=1.02) ax2.set_aspect(aspect_ratio) # Apply the corrected aspect ratio # 6. Add a common vertical colorbar on the RIGHT side # By passing ax=[ax1, ax2], Matplotlib centers it relative to both plots cbar = fig.colorbar(sc1, ax=[ax1, ax2], orientation='vertical', shrink=0.8, pad=0.03) cbar.set_label(cbar_label) plt.show() def dualplot(lats,lons,var1,var2,title1,title2,bartitle,vmin,vmax,extent,dotsize=8,date=''): center_lat = np.nanmean(lats) aspect_ratio = 1.0 / np.cos(np.radians(abs(center_lat))) # 3. Create the Map Plot fig=plt.figure(figsize=(11,8)) #fig, (ax1, ax2) = plt.subplots( 1, 2, figsize=(10, 8), # subplot_kw = { 'projection':ccrs.PlateCarree() }) # layout='constrained') gs=fig.add_gridspec(1,3,width_ratios=[4,4,0.3]) ax1=fig.add_subplot(gs[0,0], projection=ccrs.PlateCarree() ) ax2=fig.add_subplot(gs[0,1], projection=ccrs.PlateCarree() ) cax=fig.add_subplot(gs[0,2]) #ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) # Create a scatter plot of the data # We use c=statuses to color-code by the status variable #print('##',len(lons),len(var) ) scatter1 = ax1.scatter(lons, lats, c=var1, cmap='viridis', transform=ccrs.PlateCarree(), s=dotsize, edgecolor='none', zorder=5, vmin=vmin,vmax=vmax) ax1.set_title(title1+'\n'+date,color='black',y=1.02) scatter2 = ax2.scatter(lons, lats, c=var2, cmap='viridis', transform=ccrs.PlateCarree(), s=dotsize, edgecolor='none', zorder=5, vmin=vmin,vmax=vmax) ax2.set_title(title2+'\n'+date,color='black',y=1.02) for ax in [ax1, ax2]: ax.set_extent(extent) ax.set_aspect(aspect_ratio) # Apply the corrected aspect ratio ax.xaxis.set_major_locator(ticker.MultipleLocator(1)) ax.yaxis.set_major_locator(ticker.MultipleLocator(1)) ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%d°')) ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%d°')) ax.add_feature(cfeature.LAND, facecolor='lightgray', alpha=0.5) ax.add_feature(cfeature.OCEAN, facecolor='lightblue', alpha=0.5) ax.add_feature(cfeature.COASTLINE, linewidth=1) # Add a colorbar to interpret the status colors #plt.colorbar(scatter, label=bartitle, shrink=0.7) #cbar=fig.colorbar(scatter2, ax=[ax1, ax2], orientation='vertical',shrink=0.8, pad=0.1) cbar=fig.colorbar(scatter1, cax=cax, orientation='vertical',shrink=0.5, pad=0.1) # Set the map extent to frame the data properly with a little padding # Looking at the coordinates, this data is located off the coast of Southern California #padding = 0.05 #ax.set_extent([ min(lons) - padding, max(lons) + padding, # min(lats) - padding, max(lats) + padding ], crs=ccrs.PlateCarree()) # Add gridlines and coordinate labels gl1 = ax1.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False, linestyle='--', color='gray', alpha=0.7) gl1.top_labels = False gl1.right_labels = False gl2 = ax2.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False, linestyle='--', color='gray', alpha=0.7) gl2.top_labels = False gl2.right_labels = False gl2.left_labels = False plt.ion() plt.show() return(fig) def triplot(lats,lons,var1,var2,var3,title1,title2,title3,bartitle,vmin,vmax,extent,dotsize=8,date=''): center_lat = np.nanmean(lats) aspect_ratio = 1.0 / np.cos(np.radians(abs(center_lat))) # 3. Create the Map Plot fig=plt.figure(figsize=(15,8)) gs=fig.add_gridspec(1,4,width_ratios=[4,4,4,0.3]) ax1=fig.add_subplot(gs[0,0], projection=ccrs.PlateCarree() ) ax2=fig.add_subplot(gs[0,1], projection=ccrs.PlateCarree() ) ax3=fig.add_subplot(gs[0,2], projection=ccrs.PlateCarree() ) cax=fig.add_subplot(gs[0,3]) # Create a scatter plot of the data # We use c=statuses to color-code by the status variable #print('##',len(lons),len(var) ) scatter1 = ax1.scatter(lons, lats, c=var1, cmap='viridis', transform=ccrs.PlateCarree(), s=dotsize, edgecolor='none', zorder=5, vmin=vmin,vmax=vmax) ax1.set_title(title1+'\n'+date,color='black',y=1.02) scatter2 = ax2.scatter(lons, lats, c=var2, cmap='viridis', transform=ccrs.PlateCarree(), s=dotsize, edgecolor='none', zorder=5, vmin=vmin,vmax=vmax) ax2.set_title(title2+'\n'+date,color='black',y=1.02) scatter3 = ax3.scatter(lons, lats, c=var3, cmap='viridis', transform=ccrs.PlateCarree(), s=dotsize, edgecolor='none', zorder=5, vmin=vmin,vmax=vmax) ax3.set_title(title3+'\n'+date,color='black',y=1.02) for ax in [ax1, ax2, ax3]: ax.set_extent(extent) ax.set_aspect(aspect_ratio) # Apply the corrected aspect ratio ax.xaxis.set_major_locator(ticker.MultipleLocator(1)) ax.yaxis.set_major_locator(ticker.MultipleLocator(1)) ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%d°')) ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%d°')) ax.add_feature(cfeature.LAND, facecolor='lightgray', alpha=0.5) ax.add_feature(cfeature.OCEAN, facecolor='lightblue', alpha=0.5) ax.add_feature(cfeature.COASTLINE, linewidth=1) # Add a colorbar to interpret the status colors #plt.colorbar(scatter, label=bartitle, shrink=0.7) #cbar=fig.colorbar(scatter2, ax=[ax1, ax2], orientation='vertical',shrink=0.8, pad=0.1) cbar=fig.colorbar(scatter1, cax=cax, orientation='vertical',shrink=0.5, pad=0.1) # Set the map extent to frame the data properly with a little padding # Looking at the coordinates, this data is located off the coast of Southern California #padding = 0.05 #ax.set_extent([ min(lons) - padding, max(lons) + padding, # min(lats) - padding, max(lats) + padding ], crs=ccrs.PlateCarree()) # Add gridlines and coordinate labels gl1 = ax1.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False, linestyle='--', color='gray', alpha=0.7) gl1.top_labels = False gl1.right_labels = False gl2 = ax2.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False, linestyle='--', color='gray', alpha=0.7) gl2.top_labels = False gl2.right_labels = False gl2.left_labels = False gl3 = ax3.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False, linestyle='--', color='gray', alpha=0.7) gl3.top_labels = False gl3.right_labels = False gl3.left_labels = False plt.ion() plt.show() return(fig) #-----main-------------# #Read logs from 1dvar from an AMPR field experiment #Make maps of statistics like retrieval status, TB error #Find wx maps to overlay? Use Worldview? #datelabel='20241107' #datelabel='20241113' datelabel='20241018' #penfile='out/20241023_new_save/penalty_steps.txt' #penfile='out/20241023_new/penalty_steps.txt' penfile='out/'+datelabel+'/penalty_steps.txt' #status, lat, lon 3 parameters #tbpen steps #xpen steps data_array=[] init_tbpens=[] final_tbpens=[] print('open') n=-1 with open(penfile, 'r') as file: for line in file: line=line.strip() print(line) if not line: continue #skips empty lines if line.startswith('TBpen:'): tbpens=[float(x) for x in line.split()[1:]] #print('tbpens:',tbpens) #if(len(tbpens)>0): # print('tbpens:',tbpens) # init_tbpens.append(tbpens[0]) # final_tbpens.append(tbpens[-1]) #else: # init_tbpens.append(-999) # final_tbpens.append(-999) elif line.startswith('Xpen:'): xpens=[float(x) for x in line.split()[1:]] #print('xpens):',xpens) #if(len(xpens)>0): # init_xpens.append(xpens[0]) # final_xpens.append(xpens[-1]) #else: # init_xpens.append(0) # final_xpens.append(999.) #xpen is last line for this point, so define the dict here. elif line.startswith('tbinit:'): tbinit=[float(x) for x in line.split()[1:]] #print('tbinit):',tbinit) elif line.startswith('tbvec:'): tbvec=[float(x) for x in line.split()[1:]] #print('tbvec):',tbvec) elif line.startswith('tbmerra:'): tbmerra=[float(x) for x in line.split()[1:]] #print('tbmerra):',tbmerra) tbdiffs=np.array(tbinit)-np.array(tbvec) #tbpen_recalc=np.sqrt(np.sum(tbdiffs)**2/2.25/len(tbdiffs)) #old one tbpen_recalc=np.sqrt(np.sum(tbdiffs**2/2.25)/len(tbdiffs)) row_dict={ 'lat': lat, 'lon': lon, 'pos': pos, 'eia': eia, 'status': status, 'iwv': iwv, 'clw': clw, 'sws': sws, 'tbpens': tbpens, #a list 'xpens': xpens, #a list 'tbpen_recalc': tbpen_recalc, 'n_iter': len(tbpens), 'm_tiwv': m_tiwv, 'm_clw': m_clw, 'm_sws': m_sws, 'm_tbpen': m_tbpen, 'tbinit': tbinit, 'tbvec': tbvec, 'tbmerra': tbmerra, } #if n>2000: # print('exiting early for testing') # continue data_array.append(row_dict) elif line.startswith('Merra:'): print('Merra') parts=[float(x) for x in line.split()[1:]] m_tiwv=parts[0] m_clw=parts[1] m_sws=parts[2] m_tbpen=parts[3] else: n=n+1 print('other') parts = line.split() #print(parts) # Ensure the line has the expected columns before parsing if len(parts) >= 6: status = float(parts[0]) lat = float(parts[1]) lon = float(parts[2]) eia = float(parts[3]) pos = float(parts[4]) iwv = float(parts[5]) #mm clw = float(parts[6]) #m/s? sws = float(parts[7]) #m/s? print('after loop') #check sizes #print() #vars_size = [ # (name, sys.getsizeof(value)) # for name, value in locals().items() # if not name.startswith('_') and name not in ['sys', 'vars_size'] #] # ## Sort by size (largest first) and print #for name, size in sorted(vars_size, key=lambda x: x[1], reverse=True): # print(f"{name:<15} : {size} bytes") #print() #truncate if datasets are unequal (file still in progress?): l1=len(data_array)-1 #Avoid last one in case incomplete #l2=len(final_tbpens) #print('l1,l2',l1,l2) #print('setting arrays to length ',l2) #lats=lats[0:l2] #lons=lons[0:l2] #final_tbpens=final_tbpens[0:l2] #statuses=statuses[0:l2] #eias=eias[0:l2] #swss=swss[0:l2] #iwvss=iwvs[0:l2] #print('lens',len(lons),len(lats),' ',len(final_tbpens),len(statuses),len(eias)) data_array=data_array[0:l1] lats=[ row['lat'] for row in data_array ] lons=[ row['lon'] for row in data_array ] init_tbpens=[ row['tbpens'][0] for row in data_array ] init_xpens=[ row['xpens'][0] for row in data_array ] final_tbpens=[ row['tbpens'][-1] for row in data_array ] final_xpens=[ row['xpens'][-1] for row in data_array ] statuses=np.array([ row['status'] for row in data_array ]) #must use np.array to take subset with _where_ eias=[ row['eia'] for row in data_array ] iwvs=np.array([ row['iwv'] for row in data_array ]) clws=np.array([ row['clw'] for row in data_array ]) swss=np.array([ row['sws'] for row in data_array ]) tbpen_recalc=[ row['tbpen_recalc'] for row in data_array ] m_tbpens=[ row['m_tbpen'] for row in data_array ] m_iwvs=[ row['m_tiwv'] for row in data_array ] m_swss=[ row['m_sws'] for row in data_array ] m_clws=[ row['m_clw'] for row in data_array ] tbinits=np.array([ row['tbinit'] for row in data_array ]) tbvecs=np.array([ row['tbvec'] for row in data_array ]) tbmerras=np.array([ row['tbmerra'] for row in data_array ]) print(len(final_tbpens)) #extent=[-119.2,-118,26,28] #extent=[-118.7,-117.8,25.25,26.25] #extent=[-120,-117,22.5,28] extent=[min(lons)-.5,max(lons)+.5,min(lats)-.5,max(lats)+.5] ww=np.where(statuses<0)[0] #title='Retrieval TB Penalty' #bartitle='TBpen' #vmin=0 #vmax=10 #myplot(lats,lons,final_tbpens,title,bartitle,vmin,vmax,extent,date=datelabel) #title='tbinit ch1' #bartitle='TB [K]' #vmin=115 #vmax=150 #tb1init=[tbinits[k][0] for k in range(0,len(tbinits))] #myplot(lats,lons,tb1init,title,bartitle,vmin,vmax,extent,date=datelabel) # #title='tbinit ch2' #bartitle='TB [K]' #vmin=110 #vmax=125 #tb2init=[tbinits[k][1] for k in range(0,len(tbinits))] #myplot(lats,lons,tb2init,title,bartitle,vmin,vmax,extent,date=datelabel) # #title='tbinit ch3' #bartitle='TB [K]' #vmin=130 #vmax=160 #tb3init=[tbinits[k][2] for k in range(0,len(tbinits))] #myplot(lats,lons,tb3init,title,bartitle,vmin,vmax,extent,date=datelabel) title1='AMPR Tb 10V' title2='AMPR Tb 10H' bartitle='TB [K]' vmin=90 vmax=150 dualplot(lats,lons,tbvecs[:,0],tbvecs[:,1],title1,title2,bartitle, vmin,vmax,extent,date=datelabel) title1='AMPR Tb 19V' title2='AMPR Tb 19H' bartitle='TB [K]' vmin=100 vmax=150 dualplot(lats,lons,tbvecs[:,2],tbvecs[:,3],title1,title2,bartitle, vmin,vmax,extent,date=datelabel) title1='AMPR Tb 37V' title2='AMPR Tb 37H' bartitle='TB [K]' vmin=120 vmax=210 dualplot(lats,lons,tbvecs[:,4],tbvecs[:,5],title1,title2,bartitle, vmin,vmax,extent,date=datelabel) title1='AMPR Tb 85V' title2='AMPR Tb 85H' bartitle='TB [K]' vmin=160 vmax=220 dualplot(lats,lons,tbvecs[:,6],tbvecs[:,7],title1,title2,bartitle, vmin,vmax,extent,date=datelabel) #title='tbmerra ch1' #bartitle='TB [K]' #vmin=115 #vmax=150 #tb1merra=[tbmerras[k][0] for k in range(0,len(tbinits))] #myplot(lats,lons,tb1merra,title,bartitle,vmin,vmax,extent,date=datelabel) title='Retrieval Status' bartitle='Status (>=0 is good)' vmin=-5 vmax=5 print(len(statuses)) myplot(lats,lons,statuses,title,bartitle,vmin,vmax,extent,date=datelabel) ## title='Incidence Angles' ## bartitle='Incidence Angle (Degrees)' ## vmin=0 ## vmax=70 ## myplot(lats,lons,eias,title,bartitle,vmin,vmax,extent) title1='Integrated Water Vapor' title2='Merra Integrated Water Vapor' bartitle='IWV (mm)' vmin=0 vmax=20 iwvs[ww]=np.nan #myplot(lats,lons,iwvs,title1,bartitle,vmin,vmax,extent,date=datelabel) bartitle='IWV (mm)' #myplot(lats,lons,m_iwvs,title2,bartitle,vmin,vmax,extent,date=datelabel) dualplot(lats,lons,iwvs,m_iwvs,title1,title2,bartitle, vmin,vmax,extent,date=datelabel) title1='Surface Wind Speed' title2='Merra Surface Wind Speed' bartitle='Wind Speed (m/s)' vmin=0 vmax=15 swss[ww]=np.nan dualplot(lats,lons,swss,m_swss,title1,title2,bartitle, vmin,vmax,extent,date=datelabel) title1='AMPR Retrieved Cloud Liquid Water' title2='Merra Cloud Liquid Water' bartitle='CLW (mg/m3)' vmin=0 vmax=30 clws[ww]=np.nan dualplot(lats,lons,clws,m_clws,title1,title2,bartitle, vmin,vmax,extent,date=datelabel) title1='Initial TBpen' title2='Merra TBpen' bartitle='Avg Norm TB error' vmin=0 vmax=10 dualplot(lats,lons,init_tbpens,m_tbpens,title1,title2,bartitle, vmin,vmax,extent,date=datelabel) triplot(lats,lons, init_tbpens,final_tbpens,m_tbpens, 'Initial TBpen','Final TBpen','Merra TBpen',bartitle, vmin,vmax,extent,date=datelabel) #title='Initial TBpen Recalc' #bartitle='Avg Norm TB error' #vmin=0 #vmax=10 #myplot(lats,lons,tbpen_recalc,title,bartitle,vmin,vmax,extent,date=datelabel) #make it log? #CBB also plot Merra SWS, IWV #add these to the retrieval code plt.tight_layout() pdb.set_trace()