#!/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 pickle 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)) fig = plt.figure(figsize=(8,6)) # 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.4, pad=0.1, label=bartitle) # 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? savefile='out/20241104/savefile.20241104.pkl' datelabel='20241104' #savefile='out/20241113/savefile.20241113.pkl' #datelabel='20241113' #savefile='out/20241107/savefile.20241107.pkl' #does not exist yet #datelabel='20241107' #logfile='out/20241023_new_save/log.txt' # #lat lon status # #OR lat lon status, 3 parameters # ##penfile='out/20241023_new_save/penalty_steps.txt' #always 3 lines #penfile='out/20241023_new/penalty_steps.txt' #always 3 lines # #status, lat, lon 3 parameters # #tbpen steps # #xpen steps # #stepfile='out/20241023_new_save/steps.txt' # #x,y,eia # #possibly: # #initial xarr # #inital xmerra # #xarr history # #lat, lon, status data_array=[] print('open') with open(savefile,'rb') as file: fd=pickle.load(file) #flight data #extent=[-120,-117,22.5,28] #extent=[-119.2,-118,26,28] #extent=[-118.7,-117.8,25.25,26.25] lats=np.array(fd['lat']) lons=np.array(fd['lon']) ww=np.where(lats==0.0) lats[ww]=np.nan lons[ww]=np.nan tbobs=np.array(fd['tbobs']) #(50,4411,8) tbmerra=np.array(fd['tbmerra']) statuses=np.array(fd['status']) iwvs=np.array(fd['iwv']) m_iwvs=np.array(fd['m_iwv']) swss=np.array(fd['sws']) m_swss=np.array(fd['m_sws']) clws=np.array(fd['clw']) m_clws=np.array(fd['m_clw']) iwvs[ww]=np.nan swss[ww]=np.nan clws[ww]=np.nan uu=np.where(statuses<0) iwvs[uu]=np.nan swss[uu]=np.nan clws[uu]=np.nan #Need to deal with default value of zeros for unfinished files. extent=[np.nanmin(lons.flatten()),np.nanmax(lons.flatten()), np.nanmin(lats.flatten()),np.nanmax(lats.flatten())] # title='Retrieval TB Penalty' # bartitle='TBpen' # vmin=0 # vmax=10 # myplot(lats,lons,fd['tbpen'],title,bartitle,vmin,vmax,extent) # # title='tbinit ch1' # bartitle='TB [K]' # vmin=115 # vmax=150 # #tb1init=[tbinits[k][0] for k in range(0,len(tbinits))] # myplot(lats,lons,fd['tbinit'][:,:,0],title,bartitle,vmin,vmax,extent) # # title='tbinit ch2' # bartitle='TB [K]' # vmin=110 # vmax=125 # #tb2init=[tbinits[k][1] for k in range(0,len(tbinits))] # myplot(lats,lons,fd['tbinit'][:,:,1],title,bartitle,vmin,vmax,extent) # # title='tbinit ch3' # bartitle='TB [K]' # vmin=130 # vmax=160 # #tb3init=[tbinits[k][2] for k in range(0,len(tbinits))] # myplot(lats,lons,fd['tbinit'][:,:,2],title,bartitle,vmin,vmax,extent) # # title='tbobs ch1' # bartitle='TB [K]' # vmin=115 # vmax=150 # #tb1vec=[tbvecs[k][0] for k in range(0,len(tbinits))] # myplot(lats,lons,fd['tbobs'][:,:,0],title,bartitle,vmin,vmax,extent) # # title='tbmerra ch1' # bartitle='TB [K]' # vmin=115 # vmax=150 # #tb1merra=[tbmerras[k][0] for k in range(0,len(tbinits))] # myplot(lats,lons,fd['tbmerra'][:,:,0],title,bartitle,vmin,vmax,extent) # # title='Retrieval Status' # bartitle='Status (>=0 is good)' # vmin=-5 # vmax=5 # #print(len(statuses)) # myplot(lats,lons,fd['status'],title,bartitle,vmin,vmax,extent) # # ## title='Incidence Angles' # ## bartitle='Incidence Angle (Degrees)' # ## vmin=0 # ## vmax=70 # ## myplot(lats,lons,eias,title,bartitle,vmin,vmax,extent) # # title='Integrated Water Vapor' # bartitle='IWV (mm)' # vmin=0 # vmax=20 # myplot(lats,lons,fd['iwv'],title,bartitle,vmin,vmax,extent) # # title='Merra Integrated Water Vapor' # bartitle='IWV (mm)' # vmin=0 # vmax=30 # myplot(lats,lons,fd['m_iwv'],title,bartitle,vmin,vmax,extent) # # title='Surface Wind Speed' # bartitle='Wind Speed (m/s)' # vmin=0 # vmax=10 # myplot(lats,lons,fd['sws'],title,bartitle,vmin,vmax,extent) # # title='Merra Surface Wind Speed' # bartitle='Wind Speed (m/s)' # vmin=0 # vmax=10 # myplot(lats,lons,fd['m_sws'],title,bartitle,vmin,vmax,extent) # # title='Cloud Liquid Water' # bartitle='CLW (mg/m3)' # vmin=0 # vmax=10 # myplot(lats,lons,fd['clw'],title,bartitle,vmin,vmax,extent) # # title='Merra Cloud Liquid Water' # bartitle='CLW (mg/m3)' # vmin=0 # vmax=10 # myplot(lats,lons,fd['m_clw'],title,bartitle,vmin,vmax,extent) # # # title='Merra TBpen' # bartitle='Avg Norm TB error' # vmin=0 # vmax=5 # myplot(lats,lons,fd['mtbpen'],title,bartitle,vmin,vmax,extent) #make it log? # # #will need to calculate # #title='Initial TBpen' # #bartitle='Avg Norm TB error' # #vmin=0 # #vmax=5 # #myplot(lats,lons,init_tbpen,title,bartitle,vmin,vmax,extent) #make it log? # # # #CBB also plot Merra SWS, IWV # #add these to the retrieval code # # # # pdb.set_trace() # print() print() title1='AMPR Tb 10V' title2='AMPR Tb 10H' bartitle='TB [K]' vmin=90 vmax=180 lats=lats.flatten() lons=lons.flatten() tb0=tbobs[:,:,0].flatten() tb1=tbobs[:,:,1].flatten() tb2=tbobs[:,:,2].flatten() tb3=tbobs[:,:,3].flatten() tb4=tbobs[:,:,4].flatten() tb5=tbobs[:,:,5].flatten() tb6=tbobs[:,:,6].flatten() tb7=tbobs[:,:,7].flatten() dualplot(lats,lons,tb0,tb1,title1,title2,bartitle, vmin,vmax,extent,date=datelabel) title1='AMPR Tb 19V' title2='AMPR Tb 19H' bartitle='TB [K]' vmin=90 vmax=130 dualplot(lats,lons,tb2,tb3,title1,title2,bartitle, vmin,vmax,extent,date=datelabel) title1='AMPR Tb 37V' title2='AMPR Tb 37H' bartitle='TB [K]' vmin=90 vmax=140 dualplot(lats,lons,tb4,tb5,title1,title2,bartitle, vmin,vmax,extent,date=datelabel) title1='AMPR Tb 85V' title2='AMPR Tb 85H' bartitle='TB [K]' vmin=90 vmax=150 dualplot(lats,lons,tb6,tb7,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=1 print(len(statuses)) myplot(lats,lons,statuses.flatten(),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=30 #myplot(lats,lons,iwvs.flatten(),title1,bartitle,vmin,vmax,extent,date=datelabel) bartitle='IWV (mm)' vmin=0 vmax=30 #myplot(lats,lons,m_iwvs,title2,bartitle,vmin,vmax,extent,date=datelabel) dualplot(lats,lons,iwvs.flatten(),m_iwvs.flatten(),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=20 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=80 dualplot(lats,lons,clws,m_clws,title1,title2,bartitle, vmin,vmax,extent,date=datelabel) #init_tbpen will have to be calculated for this one. #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) pdb.set_trace()