#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Thu Sep 15 08:22:11 2022 @author: badler This script creates varies type of meteorological and housekeeping plots It uses data from the MWR, ASSIST, MET station, and TROPoe It is called from master_tropoe.csh. The variable 'combination' determines which types of plots are created """ import netCDF4 as nc import os import glob import matplotlib.dates as mdates import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import make_axes_locatable from matplotlib import gridspec from cycler import cycler from scipy.interpolate import interp1d from metpy.units import units import numpy.matlib import metpy.calc as mpcalc import datetime import numpy as np import pandas as pd from pandas.tseries.frequencies import to_offset import xarray as xr from netCDF4 import Dataset import scipy.io as sio import shutil #plot settings params = {'legend.fontsize': 'xx-large', 'figure.figsize': (15, 5), 'axes.labelsize': 'x-large', 'axes.titlesize':'xx-large', 'xtick.labelsize':'x-large', 'ytick.labelsize':'x-large', 'xtick.minor.size':3.5, 'xtick.top': True, 'ytick.right': True} plt.rcParams.update(params) default_cycler=(cycler(color=plt.cm.nipy_spectral(np.linspace(0,1,15)))) plt.rc('axes', prop_cycle=default_cycler) #some functions def datenum_to_datetime(datenum): """ from: https://gist.github.com/victorkristof/b9d794fe1ed12e708b9d Convert Matlab datenum into Python datetime. :param datenum: Date in datenum format :return: Datetime object corresponding to datenum. """ days = datenum % 1 return datetime.datetime.fromordinal(int(datenum)) \ + datetime.timedelta(days=days) \ - datetime.timedelta(days=366) def datestr2num(*args): x=args[0] if len(args)==1: str_format='%Y%m%d' else: str_format=args[1] n=mdates.date2num(datetime.datetime.strptime(x,str_format)) return n def num2datestr(*args): x = args[0] if len(args)==1: str_format='%Y%m%d%H%M%S' else: str_format=args[1] d=mdates.num2date(x) s=d.strftime(str_format) return s def num2doy(datestring): date_object = datetime.datetime.strptime(datestring, "%Y%m%d") return date_object.timetuple().tm_yday def get_ticks(start, end): from datetime import timedelta as td delta = end - start if delta <= td(days=1.5): loc = mdates.HourLocator(byhour=range(0,24,6)) fmt = mdates.DateFormatter('%H:%M') elif delta <= td(days=3): loc = mdates.HourLocator(byhour=range(0,24,12)) fmt = mdates.DateFormatter('%m/%d %H') elif delta <= td(days=30): loc = mdates.DayLocator() fmt = mdates.DateFormatter('%m/%d') elif delta <= td(days=100): loc = mdates.DayLocator(interval=5) fmt = mdates.DateFormatter('%m/%d') else: loc = mdates.DayLocator(interval=20) fmt = mdates.DateFormatter('%m/%d') return loc,fmt def calc_pressureprofile(p0,temp,height): #function to calculate pressure profile from surface pressure # and temperature profiles using barometric height formulas #temperature in K, pressure in hPa, output in hPa pres=np.full(temp.shape,np.nan) pres[0] = p0 for z in range(1,len(height)): pres[z] = pres[z-1]*np.exp(-9.81*(height[z]-height[z-1])/(287 * np.nanmean([temp[z],temp[z-1]]))) return pres def calc_theta(temp,pres): #function to calculate potential tempeature difference #temperature in K, pressure in hPa theta = temp*(1000./pres)**0.2854 return theta def get_nearest(refa,mya): #function to get idx of values in mya which are closest to values in refa idx=abs(refa[:,None]-mya[None,:]).argmin(axis=-1) return idx def vaporpressure(temp): #calculates vapor pressure over a flat water surface #if the input temperature is the air temperature it returns # saturation vapor pressure #if input is dew point temperature #it returns actual vapor pressure #tempeature in celsius wv = 6.112*np.exp(17.62*temp/(243.12+temp)) return wv def calc_mixingratio(e,p): #e: in hPa, p in hpa #return mr in kg/kg mr=0.622*(e/(p-e)) return mr def calc_absolutehumidity(e,temp): #function to calculate absolute humidity # from actual water vapour e (hPa) and air temperature (C) #output in g/m^3 a = 100.*e/(461.52*(temp+273.15))*1000 return a def calc_iwv(abshum,height): #compute iwv from profile, abshum in g/m3 and height in m #compute mean of abshum for height layers #previous method, abshum and diff(hegiht) don't have the same shape! # abshummean=moving_average(abshum,2) # iwv=np.nansum(abshummean*np.diff(height))/1000 df=pd.DataFrame(np.transpose(abshum/1000)) dfmean=df.rolling(2).mean() amean=np.transpose(dfmean.to_numpy()[1:,:]) iwv=np.nansum(amean*np.diff(height),axis=1) return iwv #%% #get some global variables funcdir=os.getenv('funcdir') campaigndir=os.getenv('campaigndir') combination=os.getenv('combination') datestring=os.getenv('pdate') site=os.getenv('site') oversion=os.getenv('oversion') fpathout=campaigndir+'/quicklooks' #xmin=datestr2num(datestring) #xmax=datestr2num(datestring)+1 #24 hr periods sdate=pd.to_datetime(datestring) edate=pd.to_datetime(datestring)+pd.to_timedelta('1D') #for realtime processing #if os.environ.get('statustime'): # edate=pd.Timestamp(os.environ.get('statustime')) # sdate=edate-pd.to_timedelta('1D') #else: # sdate=pd.Timestamp.now()-pd.to_timedelta('1D') # edate=pd.Timestamp.now() xmin=mdates.date2num(sdate) xmax=mdates.date2num(edate) yesterday=num2datestr(datestr2num(datestring)-1,'%Y%m%d') plotinfo={} #some general defintions tsdefcolor=['black','indianred','forestgreen','royalblue','mediumpurple','magenta','cyan','lime','orange'] tscycler=cycler(color=tsdefcolor) #read data which are needed regardless of 'combination' #%% #MET data that was used as input listdir = sorted(glob.glob(os.path.join(campaigndir,'met','*'+datestring+'*.nc'))) if len(listdir) == 0: print('no netcdf MET data found for '+datestring) else: print('read '+listdir[0]) IN=xr.open_dataset(listdir[0]) dt=pd.to_datetime(IN.base_time.values+IN.time_offset.values,unit='s') IN['time']=dt if 'atmosphericPressure' in IN.keys(): IN['atmos_pressure']=IN.atmosphericPressure elif 'pressure' in IN.keys(): IN['atmos_pressure']=IN.pressure if 'temperature' in IN.keys(): IN['temp_mean']=IN['temperature'] if 'relativeHumidity' in IN.keys(): IN['rh_mean']=IN['relativeHumidity'] elif 'relative_humidity' in IN.keys(): IN['rh_mean']=IN['relative_humidity'] #IN['atmos_pressure']=IN['atmos_pressure']*10 #is in kPa #compute water vapor mixing ratio E=vaporpressure(IN.temp_mean) e=IN.rh_mean/100*E IN['mr_mean']=calc_mixingratio(e,IN.atmos_pressure) #convert to g/kg IN['mr_mean']=IN['mr_mean']*1000 metpsl=IN.copy() ##%% #ceilometer data listdir = sorted(glob.glob(os.path.join(campaigndir,'cbh','*'+datestring+'*.cdf'))) if len(listdir) == 0: print('no ceilometer data found for '+datestring) else: print('read '+listdir[0]) IN=xr.open_dataset(listdir[0]) dt=pd.to_datetime(IN.base_time.values+IN.time_offset.values,unit='s') IN['time']=dt try: IN['cbh']=IN['cloudHeight'] except: IN['cbh']=IN['first_cbh'] try: if IN['cbh'].units == 'm': unit='m' udiv=1000 except: print('no info on ceilo cbh unit,assume it is km') udiv=1 #make sure cbh is float IN['cbh']=IN.cbh.astype(float)/udiv IN['cbh'][IN['cbh']<0]=np.nan # IN.cloudHeight=IN.cloudHeight*1000 ceil=IN.copy() ##%% #read present day listdir = sorted(glob.glob(os.path.join(campaigndir,'retrieval_output','*tropoe'+combination+'*.'+oversion+'*.*'+datestring+'*.nc'))) try: listdiry = sorted(glob.glob(os.path.join(campaigndir,'retrieval_output','*tropoe'+combination+'*.'+oversion+'*.*'+yesterday+'*.nc'))) listdir=listdiry+listdir except: print('no TROPoe data for yesterday') if len(listdir)==0: print(os.path.join(campaigndir,'retrieval_output','*tropoe'+combination+'*.'+oversion+'*.*'+datestring+'*.nc')) print('no TROPoe output found for version '+oversion+' on '+datestring+'--> try any other version') listdir = sorted(glob.glob(os.path.join(campaigndir,'retrieval_output','*tropoe*.'+datestring+'*.nc'))) if len(listdir)>1: print('now I found more than on version, I do not know what to do ') print(listdir) elif len(listdir)==0: print('still no TROPoe output found' ) else: #some variables have to be dropped because xarray cannot handle dimensions if 'IN' in locals(): del IN for l in listdir: print('read '+l) IN_=xr.load_dataset(l,decode_times=False) dt=pd.to_datetime(IN_.base_time.values+IN_.time_offset.values,unit='s') IN_['time']=dt if not 'IN' in locals(): IN=IN_ else: IN=xr.concat([IN,IN_], dim='time') #make it timezone aware (not sure) #IN['time']=pd.to_datetime(IN['time'].values) #IN=xr.open_dataset(listdir[0],decode_times=False,drop_variables=['dfs','Sa','Sop','Akernal']) #dt=[datetime.datetime.fromtimestamp(t,datetime.UTC) for t in IN['base_time'].values+IN['time_offset'].values] #IN['time']=dt #time stamps are not always sorted! IN=IN.sortby('time') #add index_dim as coordinate IN['index_dim']=['pwv','pblh','sbih','sbim','sbLCL'] ##%% #compute CAPE and CIN using metpy t0=IN.temperature.sel(height=0).values * units.degC t=IN.temperature.values * units.degC p=IN.pressure.values * units.mbar td0=IN.dewpt.sel(height=0).values * units.degC td=IN.dewpt.values * units.degC cape=np.full(len(t0),np.nan) cin=np.full(len(t0),np.nan) for i_t in range(len(t0)): parcel_profile=mpcalc.parcel_profile(p[i_t,:],t0[i_t],td0[i_t]).to('degC') cape_,cin_=mpcalc.cape_cin(p[i_t,:],t[i_t,:],td[i_t,:],parcel_profile,which_lfc='bottom',which_el='top') cape[i_t]=cape_.magnitude cin[i_t]=cin_.magnitude da=xr.DataArray(data=cape,dims='time',coords=dict(time=(['time'],IN.time.values))) IN['cape']=da da=xr.DataArray(data=-1*100*cin,dims='time',coords=dict(time=(['time'],IN.time.values))) IN['cin100']=da #flag everything with gamma > 3 mask=IN.gamma>3 mask=mask.compute() for f in IN.keys(): if 'time' in IN[f].dims and IN[f].dtype == np.float32 and (not 'amma' in f and not 'rms' in f and not 'base_time' in f \ and not 'time' in f): IN[f][mask]=np.nan #resample to tres, so that gaps are plotted as nan tres=float(IN.attrs['vip_tres'].split(' ')[0]) IN=IN.resample(time=str(round(tres))+'min').nearest(tolerance='1min') ##%% tropoe=IN.copy() #%% # #now read data which are combination specific if 'MP' in combination: #set tropoe to nan, where rmsr = 0 meaning that no mwr is available if not 'assist' in combination: mask = tropoe.rmsr != 0 tropoe = tropoe.where(mask) #MWR lv1 files try: mwrname=combination[combination.index('MP3'):combination.index('MP3')+7] except: mwrname=combination[combination.index('MP'):combination.index('MP')+7] #read all days #should be lv1!!! listdir = sorted(glob.glob(os.path.join(campaigndir,mwrname,'lv1',siteid+'*'+mwrname+'*.cdf'))) if len(listdir) ==0: print('no MWR data found for '+datestring) else: #and limit to all days before processed day ldates=[pd.to_datetime(l[l.index('4TROPoe')+8:l.index('4TROPoe')+16]) for l in listdir] listdir=[listdir[i] for i in range(len(ldates)) if ldates[i]=edate-datetime.timedelta(days=90)] if 'OUT' in locals(): del OUT for l in listdir: print('read '+l) IN=xr.open_dataset(l) #zenith only zenithmask=abs(IN.elev-90)<1e-5 IN=IN.where(zenithmask.compute(),drop=True) IN['t_sfc']=IN['t_sfc']-273.15 #compute water vapor mixing ratio E=vaporpressure(IN.t_sfc) e=IN.rh_sfc/100*E IN['mr_sfc']=calc_mixingratio(e,IN.p_sfc) #convert to g/kg IN['mr_sfc']=IN['mr_sfc']*1000 if not 'OUT' in locals(): OUT=IN else: OUT=xr.concat([OUT,IN], dim='time') mwr=OUT.copy() #read tb obs and for all days before present day listdir = sorted(glob.glob(os.path.join(tropoedir,'retrieval_output','*tropoe'+combination+'.'+oversion+'.*.nc'))) ldates=[] for l in listdir: if not '.v' in l: continue l_=l[l.index('.v')+1:] ld=l_[l_.index('.')+1:l_.index('.')+9] ldates.append(pd.to_datetime(ld)) #version dependendt try: ldates=[pd.to_datetime(l[l.index(oversion+'.')+len(oversion)+1:l.index(oversion+'.')+len(oversion)+9]) for l in listdir if oversion in l] except: ldates=[pd.to_datetime(l[l.index(oversion+'.')+len(oversion)+4:l.index(oversion+'.')+len(oversion)+12]) for l in listdir if oversion in l] if not oversion == 'rt': ldates_=[pd.to_datetime(l[l.index('rt'+'.')+len('rt')+1:l.index('rt'+'.')+len('rt')+9]) for l in listdir if 'rt.' in l] ldates=ldates_+ldates listdir=[listdir[i] for i in range(len(ldates)) if ldates[i]=edate-datetime.timedelta(days=90)] if 'OUT' in locals(): del OUT for l in listdir: print('read '+l) IN=xr.open_dataset(l,decode_times=False,drop_variables=['dfs','Sa','Sop','Akernal','arb']) dt=pd.to_datetime(IN.base_time.values+IN.time_offset.values,unit='s') #print(IN.base_time.values) #dt=[datetime.datetime.fromtimestamp(t,datetime.UTC) for t in IN['base_time'].values+IN['time_offset'].values] IN['time']=dt idx=np.where(IN.obs_flag==2)[0] if len(idx) ==0: print('no MWR data used') continue IN=IN.isel(obs_dim=idx) if not 'OUT' in locals(): OUT=IN else: OUT=xr.concat([OUT,IN], dim='time') #drop duplicates _,idx=np.unique(OUT['time'],return_index=True) OUT=OUT.isel(time=idx) #compute iwv, not necessary because alrady computed in TROPoe # e=vaporpressure(OUT.dewpt) # a=calc_absolutehumidity(e,OUT.temperature) # iwv=calc_iwv(a.values,OUT.height.values*1000) # da=xr.DataArray(data=iwv,dims='time',coords=dict(time=(['time'],OUT.time))) OUT['iwv']=OUT.pwv*10 #reindex so that gaps are filled with nan sfreq=to_offset(pd.to_timedelta(OUT.time.diff(dim='time').median().values)).freqstr dtnew=pd.date_range(start=OUT['time'][0].values,end=OUT['time'][-1].values,freq=sfreq) OUT=OUT.reindex({'time': dtnew}) tropoeobsmwr=OUT.copy(deep=True) # #%% if 'ASSIST' in combination or 'irs' in combination: #ASSIST summary files #read all days listdir = sorted(glob.glob(os.path.join(campaigndir,'irs_sum','*.cdf'))) if len(listdir) ==0: print('no ASSIST sum data found for '+datestring) else: #and limit to all days before processed day ldates=[pd.to_datetime(l[l.index('assistsummary')+14:l.index('assistsummary')+23]) for l in listdir] listdir=[listdir[i] for i in range(len(ldates)) if ldates[i]=edate-datetime.timedelta(days=2)] if 'OUT' in locals(): del OUT for l in listdir: print('read '+l) IN=xr.open_dataset(l,chunks={}) dt=pd.to_datetime(IN.base_time.values+IN.time.values,unit='ns') IN['time']=dt #time stamps are not always sorted! IN=IN.sortby('time') if not 'OUT' in locals(): OUT=IN else: OUT=xr.concat([OUT,IN], dim='time') assistsum=OUT.copy() #plots which are plotted regardless of combination if 'tropoe' in locals(): #%% #time height sections # for pctype in ['met','uncerts','dfs']: for pctype in ['met_twv','met_rhtd','met_thetathetae','dfs','uncerts','rass']: if pctype =='met': varsall=[['temperature','theta'],['waterVapor'],['rh']] cmaps=['CMRmap','terrain_r','rainbow_r'] cblabelstr=['Temperature (°C)','Mixing ratio (g kg$^{-1}$)','Relative humidity (%)'] cbextend=['both','max','neither'] clims=[[0,35],[0,20],[0,100]] fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.fig01_timeheight_tropoe_met.png') elif pctype =='met_twv': varsall=[['temperature'],['waterVapor']] cint=0.5 cmaps=['CMRmap','terrain_r'] cblabelstr=['Temperature (°C)','Mixing ratio (g kg$^{-1}$)'] cbextend=['both','max','neither'] clims=[[-5,25],[0,10]] fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.fig02_temp_wvmr.png') elif pctype =='met_rhtd': varsall=[['rh'],['dewpt']] cint=0.5 cmaps=['terrain_r','terrain_r'] cblabelstr=['Relative humidity (%)','Dew point temp (°C)'] cbextend=['neither','both'] clims=[[0,100],[-10,30]] fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.fig05_rh_dewpt.png') elif pctype =='met_thetathetae': varsall=[['theta'],['thetae']] cint=0.5 cmaps=['CMRmap','CMRmap'] cblabelstr=['Theta (K)','Theta-E (K)'] cbextend=['both','both'] clims=[[270,320],[300,350]] fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.fig04_theta_thetae.png') elif pctype == 'uncerts': varsall=[['sigma_temperature','temperature'],['sigma_waterVapor','waterVapor']] cmaps=['gist_ncar','gist_ncar'] cblabelstr=[r'1-$\sigma$ Temperature (°C)',r'1-$\sigma$ Mixing ratio (g kg$^{-1}$)'] cbextend=['max','max'] clims=[[0,4],[0,4]] fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.figxx_timeheight_tropoe_uncerts.png') elif pctype == 'dfs': varsall=[['cdfs_temperature','temperature'],['cdfs_waterVapor','waterVapor']] cmaps=['gist_ncar','gist_ncar'] cblabelstr=['Cum. deg of freedom\n Temperature','Cum. deg of freedom\n Mixing ratio'] cbextend=['max','max'] clims=[[0,5],[0,5]] fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.figxx_timeheight_tropoe_cdf.png') elif pctype == 'rass': varsall=[['tv'],['sigmatv']] cmaps=['CMRmap','rainbow_r'] cblabelstr=['RASS Tv (°C)','1-sigma Tv (°C)'] cbextend=['both','max'] clims=[[-5,25],[0,1.5]] fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.figxx_timeheight_tropoe_rass.png') if 'time' in tropoe.obs_flag.dims: fidx = np.where(tropoe.obs_flag.isel(time=0) == 3) rheight = tropoe.obs_dimension.isel(obs_dim=fidx[0],time=0) else: fidx = np.where(tropoe.obs_flag == 3) rheight = tropoe.obs_dimension.isel(obs_dim=fidx[0]) if len(fidx[0]) == 0: print('no rass data found in observation vector') continue rass = tropoe.obs_vector.isel(obs_dim = fidx[0]).to_dataset(name='tv') rass['tv'] = xr.where(rass['tv']==-999,np.nan,rass['tv']) rass['height'] = rheight rass['sigmatv'] = tropoe.obs_vector_uncertainty.isel(obs_dim = fidx[0]) fig,axs = plt.subplots(nrows=len(varsall),ncols=1,figsize=(13.5,4.5*len(varsall)),\ #fig,axs = plt.subplots(nrows=len(varsall),ncols=1,figsize=(10.5,3.5*len(varsall)),\ facecolor='w',edgecolor='w',gridspec_kw={'hspace': 0.07,'wspace': 0.07},sharey=True,sharex=True) for i_varsall in range(0,len(varsall)): vars_=varsall[i_varsall] if type(axs) is np.ndarray: ax=axs[i_varsall] else: ax=axs ax.set_facecolor((0.9,0.9,0.9)) var=vars_[0] if not pctype == 'rass': hidx=np.where(tropoe['height']<3500/1000)[0] plotdt=mdates.date2num(tropoe.time) plotheight=1000*tropoe['height'].isel(height=hidx) plotvar=tropoe[var].isel(height=hidx).T else: plotdt = mdates.date2num(rass.time) plotheight=1000*rass.height plotvar=rass[var].T if 'ASSIST' in combination or 'irs' in combination: #everywhere ax.pcolormesh(plotdt,plotheight,plotvar,vmin=clims[i_varsall][0],\ vmax=clims[i_varsall][1],cmap=cmaps[i_varsall],alpha=0.5) #only below cloud height=np.matlib.repmat(tropoe.height,len(tropoe.time),1) height=tropoe.height.isel(height=hidx) plotvar=plotvar.values for i_t in range(len(tropoe.time)): #only when LWP if tropoe.lwp[i_t]>1: idx=np.where(height.values>tropoe.cbh.values[i_t])[0] plotvar[idx,i_t]=np.nan pc=ax.pcolormesh(plotdt,plotheight,plotvar,vmin=clims[i_varsall][0],\ vmax=clims[i_varsall][1],cmap=cmaps[i_varsall]) else: pc=ax.pcolormesh(plotdt,plotheight,plotvar,vmin=clims[i_varsall][0],\ vmax=clims[i_varsall][1],cmap=cmaps[i_varsall]) divider = make_axes_locatable(ax) cax=divider.append_axes("right", size="3%", pad=0.05) cb=plt.colorbar(pc,cax=cax,extend=cbextend[i_varsall]) cb.set_label(label=cblabelstr[i_varsall],size='x-large') #add cbh cbh=tropoe.cbh[(tropoe.lwp>1).compute()] plotdt=mdates.date2num(cbh.time) plotvar=cbh.values*1000 ax.plot(plotdt,plotvar,'ok',markerfacecolor='white') if len(vars_)==2: plotdt=mdates.date2num(tropoe.time) plotheight=1000*tropoe['height'].isel(height=hidx) plotvar=tropoe[vars_[1]].isel(height=hidx).T cmin = np.floor(np.nanmin(plotvar)) cmax = np.floor(np.nanmax(plotvar)) cint=2 levels=np.arange(cmin,cmax,cint)[:] cont=ax.contour(plotdt,plotheight,plotvar,levels,colors='black') plt.clabel(cont,levels[1::2], inline=True, fontsize=12,\ fmt='%1.0f') ax.set_ylabel('Height (m AGL)') ax.set_ylim(0,3000) ax.set_xlim(xmin,xmax) ax.set_xlabel('Time (UTC)') axs[0].set_title(datestring+', '+site+' '+combination+', TROPoe configuration:'+oversion) loc,fmt=get_ticks(mdates.num2date(xmin),mdates.num2date(xmax)) ax.xaxis.set_major_locator(loc) ax.xaxis.set_major_formatter(fmt) locminor=mdates.HourLocator() ax.xaxis.set_minor_locator(locminor) if len(varsall)==2: yofft=0.09 elif len(varsall)==3: yofft=0.12 logo=fig.add_axes([0.1,yofft,0.05,0.05],anchor='NW',zorder=1) im=plt.imread(os.path.join(funcdir,'noaa_web.png')) logo.imshow(im) logo.axis('off') plt.subplots_adjust(bottom=0.2) plt.text(0.9, yofft, 'Postprocessed retrievals created with TROPoe (Turner and Loehnert 2014, Turner and Blumberg 2019).'+\ '\nContact: Bianca.Adler@noaa.gov, Creation date: '+datetime.datetime.now(datetime.UTC).strftime('%d %b %Y %H:%M UTC'), \ fontsize=14, transform=plt.gcf().transFigure,\ ha='right') if 'ASSIST' in combination: plt.text(0.9,yofft+0.055,'Data above cloud base (circles) have to be treated with care', fontsize=10,\ transform=plt.gcf().transFigure,\ ha='right') else: plt.text(0.9,yofft+0.055,'Cloud base (circles)', fontsize=10,\ transform=plt.gcf().transFigure,\ ha='right') print('save as '+fout) plt.savefig(fout,format='png',dpi=50,bbox_inches='tight') #%% #time series TROPoe met indices and RMS for pctype in ['indices','indices_extended']: if pctype == 'indices_extended': varsall=[['rmsr','rmsa','rmsp'],['gamma'],['lwp'],['pwv'],['sbih'],['sbim'],['sbLCL','pblh'],['cape','cin100']] # varsall=[['rmsr','rmsa','rmsp'],['gamma'],['lwp'],['pwv'],['sbih'],['sbim'],['lcl','pblh','cbh']] ylabelstr=['RMS','Gamma','LWP (g m$^{-2}$)',\ 'PWV (cm)','Surface inversion \nheight (km)','Surface inversion \nstrength (°C)',\ 'Height (km)','Energy (J kg$^{-1}$)'] ylims=[[-2,10],[0,1000],[0,100],[0,7],[0,1],[0,8],[0,3],[0,4000]] fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.figxx_tropoe_indices.png') elif pctype == 'indices': #for field catalog varsall=[['lwp'],['pwv'],['sbih'],['sbim'],['pblh'],['sbLCL','cbh']] # varsall=[['rmsr','rmsa','rmsp'],['gamma'],['lwp'],['pwv'],['sbih'],['sbim'],['lcl','pblh','cbh']] ylabelstr=['LWP (g m$^{-2}$)',\ 'PWV (cm)','Surface inversion \nheight (km)','Surface inversion \nstrength (°C)',\ 'Height (km)','Height (km)'] ylims=[[0,100],[0,7],[0,1],[0,8],[0,3],[0,3]] fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.fig06_indices.png') fig,axs = plt.subplots(nrows=len(varsall),ncols=1,figsize=(13.5,3*len(varsall)),\ facecolor='w',edgecolor='w',gridspec_kw={'hspace': 0.07,'wspace': 0.00},\ tight_layout={'pad':0},sharex=True) axs=np.atleast_1d(axs) for i_varsall in range(0,len(varsall)): vars_=varsall[i_varsall] if type(axs) is np.ndarray: ax=axs[i_varsall] else: ax=axs for i_vars in range(0,len(vars_)): var=vars_[i_vars] print(var) plotdt=mdates.date2num(tropoe.time) plotvar=tropoe[var].values plotvar[plotvar==-999]=np.nan if 'cbh' in var: plotvar[tropoe.lwp<1]=np.nan ax.plot(plotdt,plotvar,'o',color=tsdefcolor[i_vars],label=var) ax.set_ylabel(ylabelstr[i_varsall]) if 'gamma' in var: ax.set_yscale('log') ax.set_ylim(ylims[i_varsall]) ax.legend(fontsize='x-large') ax.grid('on') ax.set_xlim(xmin,xmax) ax.set_xlabel('Time (UTC)') axs[0].set_title(datestring+', '+site+' '+combination+', TROPoe configuration:'+oversion) loc,fmt=get_ticks(mdates.num2date(xmin),mdates.num2date(xmax)) ax.xaxis.set_major_locator(loc) ax.xaxis.set_major_formatter(fmt) locminor=mdates.HourLocator() ax.xaxis.set_minor_locator(locminor) logo=fig.add_axes([0.1,0.04,0.05,0.05],anchor='NW',zorder=1) im=plt.imread(os.path.join(funcdir,'noaa_web.png')) logo.imshow(im) logo.axis('off') plt.text(0.9, 0.05, 'Postprocessed retrievals created with TROPoe '+\ '(Turner and Loehnert 2014, Turner and Blumberg 2019). '+\ '\nContact: Bianca.Adler@noaa.gov, Creation date: '+datetime.datetime.now(datetime.UTC).strftime('%d %b %Y %H:%M UTC'), \ fontsize=14, transform=plt.gcf().transFigure,\ ha='right') print('save as '+fout) plt.savefig(fout,format='png',dpi=50,bbox_inches='tight') #%% #profiles try: varsall=[['temperature'],['theta'],['rh'],['waterVapor']] labelstr=['Temperature (°C)','Potential temperature (K)','Relative humidity (%)', 'Water vapor (g kg$^{-1}$)'] fig,axs = plt.subplots(ncols=len(varsall),nrows=1,figsize=(5*len(varsall),5),\ facecolor='w',edgecolor='w',gridspec_kw={'hspace': 0.07,'wspace': 0.07},sharey=True) fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.figxx_profmet_tropoe.png') timeint=pd.date_range(sdate,edate,freq='2h') #find closest time stamps plottropoe=tropoe.reindex(time=timeint,method='nearest',tolerance=pd.to_timedelta('1h'),fill_value=np.nan) legendstr=plottropoe.time.dt.strftime('%H:%M UTC').values for i_varsall in range(0,len(varsall)): vars_=varsall[i_varsall] if type(axs) is np.ndarray: ax=axs[i_varsall] else: ax=axs for i_vars in range(0,len(vars_)): var=vars_[i_vars] hidx=np.where(plottropoe['height']<5500/1000)[0] p=ax.plot(plottropoe[var].isel(height=hidx).T,1000*plottropoe['height'].isel(height=hidx).T,label=legendstr) ax.grid('on') ax.set_xlabel(labelstr[i_varsall]) ax.set_ylim(0,5000) leg=fig.legend(p,legendstr,loc='lower left',ncol=7,bbox_to_anchor=(0.162,0.95,0.735,.102),\ borderaxespad=0.,handletextpad=0.1,handlelength=2,mode='expand') leg.set_title(datestring+', '+site+' '+combination+', TROPoe configuration:'+oversion,prop={'size':'xx-large'}) axs[0].set_ylabel('Height (m AGL)') logo=fig.add_axes([0.08,0.9,0.25,0.25],anchor='NW',zorder=1) im=plt.imread(os.path.join(funcdir,'noaa_web.png')) logo.imshow(im) logo.axis('off') plt.subplots_adjust(bottom=0.2) plt.text(0.9, 0., 'Postprocessed retrievals created with TROPoe (Turner and Loehnert 2014, Turner and Blumberg 2019). Contact: Bianca.Adler@noaa.gov, Creation date: '+datetime.datetime.now(datetime.UTC).strftime('%d %b %Y %H:%M UTC'), \ fontsize=14, transform=plt.gcf().transFigure,\ ha='right') print('save as '+fout) plt.savefig(fout,format='png',dpi=50,bbox_inches='tight') #%% varsall=[['sigma_temperature'],['cdfs_temperature'],['sigma_waterVapor'],['cdfs_waterVapor']] labelstr=[r'1-$\sigma$ Temperature (°C)','Cum. degree of freedom \nTemperature',\ r'1-$\sigma$ Water vapor (g kg$^{-1}$)','Cum. degree of freedom \nWater vapor'] fig,axs = plt.subplots(ncols=len(varsall),nrows=1,figsize=(5*len(varsall),5),\ facecolor='w',edgecolor='w',gridspec_kw={'hspace': 0.07,'wspace': 0.07},sharey=True) fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.figxx_profuncerts_tropoe.png') timeint=pd.date_range(sdate,edate,freq='2h') #find closest time stamps plottropoe=tropoe.reindex(time=timeint,method='nearest',tolerance=pd.to_timedelta('1h'),fill_value=np.nan) legendstr=plottropoe.time.dt.strftime('%H:%M UTC').values for i_varsall in range(0,len(varsall)): vars_=varsall[i_varsall] if type(axs) is np.ndarray: ax=axs[i_varsall] else: ax=axs for i_vars in range(0,len(vars_)): var=vars_[i_vars] hidx=np.where(plottropoe['height']<5500/1000)[0] p=ax.plot(plottropoe[var].isel(height=hidx).T,1000*plottropoe['height'].isel(height=hidx).T,label=legendstr) ax.grid('on') ax.set_xlabel(labelstr[i_varsall]) ax.set_ylim(0,5000) leg=fig.legend(p,legendstr,loc='lower left',ncol=7,bbox_to_anchor=(0.162,0.95,0.735,.102),\ borderaxespad=0.,handletextpad=0.1,handlelength=2,mode='expand') leg.set_title(datestring+', '+site+' '+combination+', TROPoe configuration:'+oversion,prop={'size':'xx-large'}) axs[0].set_ylabel('Height (m AGL)') logo=fig.add_axes([0.08,0.9,0.25,0.25],anchor='NW',zorder=1) im=plt.imread(os.path.join(funcdir,'noaa_web.png')) logo.imshow(im) logo.axis('off') plt.subplots_adjust(bottom=0.2) plt.text(0.9, 0., 'Postprocessed retrievals created with TROPoe (Turner and Loehnert 2014, Turner and Blumberg 2019). Contact: Bianca.Adler@noaa.gov, Creation date: '+datetime.datetime.now(datetime.UTC).strftime('%d %b %Y %H:%M UTC'), \ fontsize=14, transform=plt.gcf().transFigure,\ ha='right') print('save as '+fout) plt.savefig(fout,format='png',dpi=50,bbox_inches='tight') #%% except: print('error when trying to plot profiles') else: print('cannot plot profiles, because no TROPoe retrievals on day '+datestring) #%% #if MWR in combination if 'MP' in combination: #meteorological time series varsall=[['temp_meanmetpsl','temperaturetropoe','obs_vector5tropoe','t_sfcmwr'],\ ['mr_meanmetpsl','waterVaportropoe','obs_vector6tropoe','mr_sfcmwr'],\ ['rh_meanmetpsl','rhtropoe','rh_sfcmwr'],\ ['atmos_pressuremetpsl','pressuretropoe','p_sfc_origmwr'],\ ['irtmwr'],['tbstd30GHzmwr'],['lwptropoe'],['cbhceil','cbhtropoe']] ylabelstr=['Temperature (°C)','Mixing ratio (g kg$^{-1}$)',\ 'Relative humidity (%)','Pressure (mb)','Infrared temperature (K)',\ 'TB (K)','LWP (g m$^{-2}$)','Cloud-base height (km)'] ylims=[[],[],[0,100],[],[],[0,15],[0,100],[0,7]] ndays='1D' fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.figxx_mwr_tsmet.png') fig,axs = plt.subplots(nrows=len(varsall),ncols=1,figsize=(13.5,3*len(varsall)),\ facecolor='w',edgecolor='w',gridspec_kw={'hspace': 0.07,'wspace': 0.00},\ tight_layout={'pad':0},sharex=True) for i_varsall in range(0,len(varsall)): vars_=varsall[i_varsall] if type(axs) is np.ndarray: ax=axs[i_varsall] else: ax=axs for i_vars in range(0,len(vars_)): var=vars_[i_vars] labelstr='' try: if 'metpsl' in var: pc=tsdefcolor[0] labelstr='MET PSL' fvar=var[:var.index('metpsl')] plotvar=metpsl[fvar].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(metpsl['time'].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) elif 'tropoe' in var: pc=tsdefcolor[1] labelstr='TROPoe' fvar=var[:var.index('tropoe')] if 'obs_vector' in fvar: obs_flag=int(fvar[fvar.index('obs_vector')+10:]) if len(tropoe.obs_flag.dims) > 1: idx=np.where(tropoe.obs_flag.isel(time=0)==obs_flag)[0] else: idx=np.where(tropoe.obs_flag==obs_flag)[0] plotvar=tropoe['obs_vector'][:,idx].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values labelstr='obs '+labelstr pc=tsdefcolor[3] elif 'height' in tropoe[fvar].dims: plotvar=tropoe[fvar].sel(height=0,time=slice(edate-pd.to_timedelta(ndays),edate)).values else: plotvar=tropoe[fvar].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(tropoe['time'].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) elif 'mwr' in var: pc=tsdefcolor[2] labelstr=mwrname fvar=var[:var.index('mwr')] if not fvar in mwr.keys() and not var.startswith('tbstd'): continue if var.startswith('tbstd'): fidx=np.where(mwr.freq==float(var[5:var.index('GHz')]))[0] fvar='tbsky' plotvar=np.squeeze(mwr[fvar].isel(freq=fidx).rolling(time=24,center=True).std().sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) labelstr=labelstr+' std '+str(mwr.freq[fidx][0].values)+ 'GHz' elif var.startswith('tbmean'): fidx=np.where(mwr.freq==float(var[6:var.index('GHz')]))[0] fvar='tbsky' plotvar=mwr[fvar].isel(freq=fidx).rolling(time=24,center=True).mean().sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values labelstr=labelstr+' mean '+str(mwr.freq[fidx][0].values)+ 'GHz' else: plotvar=mwr[fvar].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(mwr['time'].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) elif 'ceil' in var: pc=tsdefcolor[4] labelstr=ceilometer fvar=var[:var.index('ceil')] plotvar=ceil[fvar].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(ceil['time'].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) if 'cbh' in var: marker='o' else: marker='None' ax.plot(plotdt,plotvar,label=labelstr,color=pc,marker=marker) except NameError: print('variable '+var+' cannot be plotted, probably because no data exist') ax.legend(fontsize='x-large') ax.grid('on') if len(ylims[i_varsall])>0: ax.set_ylim(ylims[i_varsall]) ax.set_ylabel(ylabelstr[i_varsall]) ax.set_xlim(xmin,xmax) ax.set_xlabel('Time (UTC)') axs[0].set_title(datestring+', '+site+' '+combination+', TROPoe configuration:'+oversion) loc,fmt=get_ticks(mdates.num2date(xmin),mdates.num2date(xmax)) ax.xaxis.set_major_locator(loc) ax.xaxis.set_major_formatter(fmt) locminor=mdates.HourLocator() ax.xaxis.set_minor_locator(locminor) logo=fig.add_axes([0.1,0.06,0.05,0.05],anchor='NW',zorder=1) im=plt.imread(os.path.join(funcdir,'noaa_web.png')) logo.imshow(im) logo.axis('off') plt.text(0.9, 0.08, 'Postprocessed retrievals created with TROPoe '+\ '(Turner and Loehnert 2014, Turner and Blumberg 2019). '+\ '\nContact: Bianca.Adler@noaa.gov, Creation date: '+datetime.datetime.now(datetime.UTC).strftime('%d %b %Y %H:%M UTC'), \ fontsize=14, transform=plt.gcf().transFigure,\ ha='right') print('save as '+fout) plt.savefig(fout,format='png',dpi=50,bbox_inches='tight') #%% #hovmoeller of TB difference forward model - obs tbdiff=tropoeobsmwr.forward_calc-tropoeobsmwr.obs_vector #last x days instantaneous values fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.figxx_mwr_tbbias_10days.png') ndays='10D' xlim=[mdates.date2num(edate-pd.to_timedelta(ndays)),mdates.date2num(edate)] cmap=plt.get_cmap('bwr') #cmin=-10.5 #cmax=10.5 cmin=-3.5 cmax=3.5 #find gaps in MWR data deltatime=mwr.time.diff(dim='time').median() tidx=np.where(mwr.time.diff(dim='time')>deltatime*4)[0] mwrgaps=xr.DataArray(data=np.full((len(tidx)),0),dims='time',coords=dict(time=(['time'],mwr.time[tidx].values))) fig, axs = plt.subplots(nrows=4, ncols=1, sharex='col', gridspec_kw={'height_ratios': [1,1,1, 3],'hspace': 0.07}, figsize=(13.5, 15.5)) #plot standard deviation at 30 GHz ax0=axs[0] fidx=np.where(mwr.freq==30)[0] plotvar=mwr['tbsky'].isel(freq=fidx).rolling(time=24,center=True).std() plotvar=plotvar.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(mwr['time'].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) ax0.plot(plotdt,plotvar) ax0.set_ylim(0,5) #plot gaps plotvar=mwrgaps.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(mwrgaps.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) ax0.plot(plotdt,plotvar,'sr') ax0.set_xlim(xlim) ax0.grid('on') ax0.set_ylabel('TB std 30 GHz (K)') divider = make_axes_locatable(ax0) cax=divider.append_axes("right", size="5%", pad=0.05) cax.axis('off') #lwp ax1=axs[1] plotvar=tropoeobsmwr.lwp.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(tropoeobsmwr.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) ax1.plot(plotdt,plotvar,'.') ax1.set_ylim(0,100) #plot gaps plotvar=mwrgaps.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(mwrgaps.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) ax1.plot(plotdt,plotvar,'sr') ax1.grid('on') ax1.set_ylabel('LWP (g m$^{-2}$)') divider = make_axes_locatable(ax1) cax=divider.append_axes("right", size="5%", pad=0.05) cax.axis('off') #iwv ax2=axs[2] plotvar=tropoeobsmwr.iwv.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(tropoeobsmwr.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) ax2.plot(plotdt,plotvar,'.') ax2.set_ylim(0,70) #plot gaps plotvar=mwrgaps.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values # plotvar=np.full(plotvar.shape,ax2.get_ylim()[0]) plotdt=mdates.date2num(mwrgaps.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) ax2.plot(plotdt,plotvar,'sr') ax2.grid('on') ax2.set_ylabel('IWV (kg m$^{-2}$)') divider = make_axes_locatable(ax2) cax=divider.append_axes("right", size="5%", pad=0.05) cax.axis('off') ax3=axs[3] ax3.set_facecolor((0.9,0.9,0.9)) plotvar=tbdiff.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(tbdiff['time'].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) obsfreq=tbdiff.obs_dim.values nanrow=np.full((1,plotvar.shape[0]),np.nan).T plotvar=np.concatenate((plotvar,nanrow),axis=1) pc=ax3.pcolormesh(plotdt,range(plotvar.shape[1]),plotvar.T,cmap=cmap,vmin=cmin,vmax=cmax) ##%% ax3.set_yticks(np.arange(plotvar.shape[1]-1)) yticklabels=[format(d,'.3f') for d in obsfreq] ax3.set_yticklabels(yticklabels,fontsize='large') ax3.set_ylim(-0.5,plotvar.shape[1]-1.5) ax3.set_ylabel('Frequency (GHz)') # plt.clim(-5,5) ax3.grid('on') vidx=np.where(obsfreq>50)[0] if len(vidx)>0: ax3.hlines(ax3.get_yticks()[vidx[0]]-0.5,ax3.get_xlim()[0],ax3.get_xlim()[1],ls='--') divider = make_axes_locatable(ax3) cax=divider.append_axes("right", size="5%", pad=0.05) cb=plt.colorbar(pc,cax=cax,extend='both') cb.set_label(label='TB forward model - TB obs (K)',size=14) loc,fmt=get_ticks(mdates.num2date(xlim[0]),mdates.num2date(xlim[-1])) ax3.xaxis.set_major_locator(loc) ax3.xaxis.set_major_formatter(fmt) locminor=mdates.HourLocator(interval=3) ax3.xaxis.set_minor_locator(locminor) ax3.set_xlabel('Time (UTC)') ax0.set_title(site+' '+combination+', TROPoe configuration:'+oversion) logo=fig.add_axes([0.05,0.06,0.05,0.05],anchor='NW',zorder=1) im=plt.imread(os.path.join(funcdir,'noaa_web.png')) logo.imshow(im) logo.axis('off') plt.text(0.9, 0.08, 'Contact: Bianca.Adler@noaa.gov, Creation date: '+datetime.datetime.now(datetime.UTC).strftime('%d %b %Y %H:%M UTC'), \ fontsize=14, transform=plt.gcf().transFigure,\ ha='right') print('save as '+fout) plt.savefig(fout,format='png',dpi=50,bbox_inches='tight') #%% # #last x days daily mean values # fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.figxx_mwr_tbbias_dailymeans.png') # # maxtimerange=mwr.time[-1].values-mwr.time[0].values # deltatime=mwr.time[-1].values.astype('datetime64[D]')-mwr.time[0].values.astype('datetime64[D]') # if pd.to_timedelta(deltatime)>pd.to_timedelta('365D'): # ndays='365D' # else: # ndays=to_offset(pd.to_timedelta(deltatime)).freqstr # # ndays='200D' # # ndays='60D' # cmap=plt.get_cmap('bwr') # cmin=-2.5 # cmax=2.5 # # #find gaps in MWR data # deltatime=mwr.time.diff(dim='time').median() # tidx=np.where(mwr.time.diff(dim='time')>deltatime*4)[0] # mwrgaps=xr.DataArray(data=np.full((len(tidx)),0),dims='time',coords=dict(time=(['time'],mwr.time[tidx].values))) # # # # #interpolate mwr data to same time axis as tropoe # mwrint=mwr.interp(time=tbdiff.time) # # #filter TB for clear sky periods # freq30=mwrint['tbsky'].isel(freq=fidx).rolling(time=24,center=True).std() # tidx=np.where(freq30<1)[0] # tbdiffclear=tbdiff.isel(time=tidx) # freq30clear=freq30.isel(time=tidx) # # #compute means for clear sky periods # tbdiffmean=tbdiffclear.resample(time='1D').mean() # freq30mean=freq30clear.resample(time='1D').mean() # #only use means if at least 20 samples available per day # avail=freq30clear.resample(time='1D').count() # tidx=np.where(avail>20)[0] # tbdiffmean=tbdiffmean.isel(time=tidx) # freq30mean=freq30mean.isel(time=tidx) # # #fill gaps with nan # dtnew=pd.date_range(start=edate-pd.to_timedelta(ndays),end=edate,freq='1D') # tbdiffmean=tbdiffmean.reindex({'time': dtnew}) # # # #compute means for LWP (regardless of clear sky) # lwpmean=tropoeobsmwr.lwp.resample(time='1D').mean() # #compute means for IWV (regardless of clear sky) # iwvmean=tropoeobsmwr.iwv.resample(time='1D').mean() # # # fig, axs = plt.subplots(nrows=4, ncols=1, sharex='col', # gridspec_kw={'height_ratios': [1,1,1, 3],'hspace': 0.07}, # figsize=(13.5, 15.5)) # #30 GHz standard deviation # ax0=axs[0] # fidx=np.where(mwrint.freq==30)[0] # plotvar=freq30mean.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values # plotdt=mdates.date2num(freq30mean.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) # # plotvar=plotvar.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values # # plotdt=mdates.date2num(mwr['time'].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) # ax0.plot(plotdt,plotvar,'.') # ax0.set_ylim(0,1) # # #plot gaps # plotvar=mwrgaps.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values # plotdt=mdates.date2num(mwrgaps.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) # ax0.plot(plotdt,plotvar,'sr') # # ax0.grid('on') # ax0.set_ylabel('TB std 30 GHz (K)') # divider = make_axes_locatable(ax0) # cax=divider.append_axes("right", size="5%", pad=0.05) # cax.axis('off') # # #lwp # ax1=axs[1] # plotvar=lwpmean.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values # plotdt=mdates.date2num(lwpmean.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) # ax1.plot(plotdt,plotvar,'.') # ax1.set_ylim(0,100) # #plot gaps # plotvar=mwrgaps.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values # plotdt=mdates.date2num(mwrgaps.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) # ax1.plot(plotdt,plotvar,'sr') # # ax1.grid('on') # ax1.set_ylabel('LWP (g m$^{-2}$)') # divider = make_axes_locatable(ax1) # cax=divider.append_axes("right", size="5%", pad=0.05) # cax.axis('off') # # #iwv # ax2=axs[2] # plotvar=iwvmean.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values # plotdt=mdates.date2num(iwvmean.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) # ax2.plot(plotdt,plotvar,'.') # ax2.set_ylim(0,70) # #plot gaps # plotvar=mwrgaps.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values # plotdt=mdates.date2num(mwrgaps.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) # ax2.plot(plotdt,plotvar,'sr') # # ax2.grid('on') # ax2.set_ylabel('IWV (kg m$^{-2}$)') # divider = make_axes_locatable(ax2) # cax=divider.append_axes("right", size="5%", pad=0.05) # cax.axis('off') # # # #tb bias plot # ax3=axs[3] # ax3.set_facecolor((0.9,0.9,0.9)) # plotvar=tbdiffmean.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values # plotdt=mdates.date2num(tbdiffmean.time.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) # nanrow=np.full((1,plotvar.shape[0]),np.nan).T # plotvar=np.concatenate((plotvar,nanrow),axis=1) # nancol=np.full((1,plotvar.shape[1]),np.nan) # plotvar=np.concatenate((plotvar,nancol),axis=0) # plotdt=np.concatenate((plotdt,np.atleast_1d(plotdt[-1])+1)) # pc=ax3.pcolormesh(plotdt,range(plotvar.shape[1]),plotvar.T,cmap=cmap,vmin=cmin,vmax=cmax) # ##%% # ax3.set_yticks(np.arange(plotvar.shape[1]-1)) # yticklabels=[format(d,'.3f') for d in obsfreq] # ax3.set_yticklabels(yticklabels,fontsize='large') # ax3.set_ylim(-0.5,plotvar.shape[1]-1.5) # ax3.set_ylabel('Frequency (GHz)') # # plt.clim(-5,5) # ax3.grid('on') # vidx=np.where(obsfreq>50)[0] # if len(vidx)>0: # ax3.hlines(ax3.get_yticks()[vidx[0]]-0.5,ax3.get_xlim()[0],ax3.get_xlim()[1],ls='--') # # divider = make_axes_locatable(ax3) # cax=divider.append_axes("right", size="5%", pad=0.05) # cb=plt.colorbar(pc,cax=cax,extend='both') # cb.set_label(label='TB forward model - TB obs (K)',size=14) # # # # xlim=[mdates.date2num(edate-pd.to_timedelta(ndays)),mdates.date2num(edate)] # ax3.set_xlim(xlim) # loc,fmt=get_ticks(mdates.num2date(xlim[0]),mdates.num2date(xlim[-1])) # ax3.xaxis.set_major_locator(loc) # ax3.xaxis.set_major_formatter(fmt) # if all(np.diff(loc())<12): # locminor=mdates.HourLocator(interval=24*1) # else: # locminor=mdates.HourLocator(interval=24*4) # ax3.xaxis.set_minor_locator(locminor) # # ax3.set_xlabel('Time (UTC)') # ax0.set_title(site+' '+combination+', TROPoe configuration:'+oversion) # # logo=fig.add_axes([0.05,0.06,0.05,0.05],anchor='NW',zorder=1) # im=plt.imread(os.path.join(campaigndir,'functions','noaa_web.png')) # logo.imshow(im) # logo.axis('off') # # plt.text(0.9, 0.08, 'Contact: Bianca.Adler@noaa.gov, Creation date: '+datetime.datetime.now(datetime.UTC).strftime('%d %b %Y %H:%M UTC'), \ # fontsize=14, transform=plt.gcf().transFigure,\ # ha='right') # # # print('save as '+fout) # plt.savefig(fout,format='png',dpi=50,bbox_inches='tight') #%% # #time series of TB # deltatime=mwr.time[-1].values.astype('datetime64[D]')-mwr.time[0].values.astype('datetime64[D]') # if pd.to_timedelta(deltatime)>pd.to_timedelta('365D'): # ndays='365D' # else: # ndays=to_offset(pd.to_timedelta(deltatime)).freqstr # # varsall=[['tbskykband'],['tbskyvband']] # fig,axs = plt.subplots(nrows=len(varsall),ncols=1,figsize=(13.5,7.25*len(varsall)),\ # facecolor='w',edgecolor='w',gridspec_kw={'hspace': 0.07,'wspace': 0.00},\ # tight_layout={'pad':0},sharex=True) # plotinfo['fig']=fig # fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.figxx_mwr_tstb.png') # # for i_varsall in range(0,len(varsall)): # vars_=varsall[i_varsall] # if type(axs) is np.ndarray: # ax=axs[i_varsall] # else: # ax=axs # for i_vars in range(0,len(vars_)): # var=vars_[i_vars] # #only vertical # #print(mwrlv1.keys()) # eidx=np.where(abs(mwr['elev']-90)<1)[0] # plotdt=mdates.date2num(mwr['time'][eidx].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) # if var =='tbskykband': # plotvar=mwr['tbsky'][eidx,0:20] # plotvar=plotvar.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values # ax.plot(plotdt,plotvar,label='MWR K-band') # ax.set_ylabel('Tb K-band (K)') # elif var =='tbskyvband': # plotvar=mwr['tbsky'][eidx,21:] # plotvar=plotvar.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values # ax.plot(plotdt,plotvar,label='MWR V-band') # ax.set_ylabel('Tb V-band (K)') # # # ax.legend(loc='center',ncol=5) # ax.grid('on') # # ax.set_xlim(xmin,xmax) # axs[-1].set_xlabel('Time (UTC)') # axs[0].set_title(datestring+', '+site+' '+combination) # # xlim=[mdates.date2num(edate-pd.to_timedelta(ndays)),mdates.date2num(edate)] # axs[0].set_xlim(xlim) # loc,fmt=get_ticks(mdates.num2date(xlim[0]),mdates.num2date(xlim[-1])) # # hours=mdates.HourLocator() # # hfmt =mdates.DateFormatter('%H%M') # ax.xaxis.set_major_locator(loc) # ax.xaxis.set_major_formatter(fmt) # locminor=mdates.HourLocator(interval=27*24) # ax.xaxis.set_minor_locator(locminor) # # logo=fig.add_axes([0.05,0.06,0.05,0.05],anchor='NW',zorder=1) # im=plt.imread(os.path.join(campaigndir,'functions','noaa_web.png')) # logo.imshow(im) # logo.axis('off') # # plt.text(0.9, 0.07, 'Contact: Bianca.Adler@noaa.gov, Creation date: '+datetime.datetime.now(datetime.UTC).strftime('%d %b %Y %H:%M UTC'), \ # fontsize=14, transform=plt.gcf().transFigure,\ # ha='right') # # # print('save as '+fout) # plt.savefig(fout,format='png',dpi=50,bbox_inches='tight') #%% if 'ASSIST' in combination or 'irs' in combination: #meteorological time series varsall=[['temp_meanmetpsl','temperaturetropoe','mean_Tb_675_680assistsum'],\ ['mr_meanmetpsl','waterVaportropoe',],\ ['rh_meanmetpsl','rhtropoe'],\ ['atmos_pressuremetpsl','pressuretropoe'],\ ['mean_Tb_985_990assistsum'],['mean_Tb_2510_2515assistsum'],\ ['hatchOpenassistsum'],\ ['lwptropoe'],['cbhtropoe','cbhceil']] # varsall=[['temp_meanmetpsl','temperaturetropoe','obs_vector5tropoe','mean_Tb_675_680assistsum']]#,\ # ['mr_meanmetpsl','waterVaportropoe','obs_vector6tropoe','mr_sfcmwr'],\ # ['rh_meanmetpsl','rhtropoe','rh_sfcmwr'],\ # ['atmos_pressuremetpsl','pressuretropoe','p_sfc_origmwr'],\ # ['irtmwr'],['tbstd30GHzmwr'],['lwptropoe'],['cbhtropoe','cbhceil']] ylabelstr=['Temperature (°C)','Mixing ratio (g kg$^{-1}$)',\ 'Relative humidity (%)','Pressure (mb)','TB 985_990 (K)','TB 2510_2515 (K)',\ 'Hatch open (1)', 'LWP (g m$^{-2}$)','Cloud-base height (km)'] ylims=[[],[],[0,100],[],[],[],[],[0,100],[],[],[]] ndays='1D' fout=os.path.join(fpathout,combination+'.'+oversion+'.'+datestring+'.figxx_assist_tsmet.png') fig,axs = plt.subplots(nrows=len(varsall),ncols=1,figsize=(13.5,3*len(varsall)),\ facecolor='w',edgecolor='w',gridspec_kw={'hspace': 0.07,'wspace': 0.00},\ tight_layout={'pad':0},sharex=True) axs=np.atleast_1d(axs) for i_varsall in range(0,len(varsall)): vars_=varsall[i_varsall] if type(axs) is np.ndarray: ax=axs[i_varsall] else: ax=axs for i_vars in range(0,len(vars_)): var=vars_[i_vars] labelstr='' #try: if 'metpsl' in var: pc=tsdefcolor[0] labelstr='MET PSL' fvar=var[:var.index('metpsl')] plotvar=metpsl[fvar].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(metpsl['time'].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) elif 'tropoe' in var: pc=tsdefcolor[1] labelstr='TROPoe' fvar=var[:var.index('tropoe')] if 'obs_vector' in fvar: obs_flag=int(fvar[fvar.index('obs_vector')+10:]) idx=np.where(tropoe.obs_flag==obs_flag)[0] plotvar=tropoe.copy(deep=True) plotvar['time']=pd.to_datetime(tropoe.time) plotvar=plotvar['obs_vector'][:,idx].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values labelstr='obs '+labelstr pc=tsdefcolor[3] elif 'height' in tropoe[fvar].dims: plotvar=tropoe.copy(deep=True) plotvar['time']=pd.to_datetime(tropoe.time) plotvar=plotvar[fvar].sel(height=0,time=slice(edate-pd.to_timedelta(ndays),edate)).values else: plotvar=tropoe.copy(deep=True) plotvar['time']=pd.to_datetime(tropoe.time) plotvar=plotvar[fvar].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values tropoe['time']=pd.to_datetime(tropoe.time) plotdt=mdates.date2num(tropoe['time'].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) elif 'assistsum' in var: pc=tsdefcolor[2] labelstr='irs' fvar=var[:var.index('assist')] if var.startswith('mean_Tb_675_680'): plotvar=assistsum[fvar]-273.15 plotvar=plotvar.sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values labelstr=labelstr+' '+fvar elif var.startswith('tbmean'): fidx=np.where(mwr.freq==float(var[6:var.index('GHz')]))[0] fvar='tbsky' plotvar=mwr[fvar].isel(freq=fidx).rolling(time=24,center=True).mean().sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values labelstr=labelstr+' mean '+str(mwr.freq[fidx][0].values)+ 'GHz' else: plotvar=assistsum[fvar].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(assistsum['time'].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) elif 'ceil' in var: pc=tsdefcolor[4] labelstr='ceilometer' fvar=var[:var.index('ceil')] plotvar=ceil[fvar].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values plotdt=mdates.date2num(ceil['time'].sel(time=slice(edate-pd.to_timedelta(ndays),edate)).values) if 'cbh' in var: marker='o' else: marker='None' ax.plot(plotdt,plotvar,label=labelstr,color=pc,marker=marker) #except NameError: # print('variable '+var+' cannot be plotted, probably because no data exist') ax.legend(fontsize='x-large') ax.grid('on') if len(ylims[i_varsall])>0: ax.set_ylim(ylims[i_varsall]) ax.set_ylabel(ylabelstr[i_varsall]) ax.set_xlim(xmin,xmax) ax.set_xlabel('Time (UTC)') axs[0].set_title(datestring+', '+site+' '+combination+', TROPoe configuration:'+oversion) loc,fmt=get_ticks(mdates.num2date(xmin),mdates.num2date(xmax)) ax.xaxis.set_major_locator(loc) ax.xaxis.set_major_formatter(fmt) locminor=mdates.HourLocator() ax.xaxis.set_minor_locator(locminor) logo=fig.add_axes([0.1,0.06,0.05,0.05],anchor='NW',zorder=1) im=plt.imread(os.path.join(funcdir,'noaa_web.png')) logo.imshow(im) logo.axis('off') plt.text(0.9, 0.08, 'Postprocessed retrievals created with TROPoe '+\ '(Turner and Loehnert 2014, Turner and Blumberg 2019). '+\ '\nContact: Bianca.Adler@noaa.gov, Creation date: '+datetime.datetime.now(datetime.UTC).strftime('%d %b %Y %H:%M UTC'), \ fontsize=14, transform=plt.gcf().transFigure,\ ha='right') print('save as '+fout) plt.savefig(fout,format='png',dpi=50,bbox_inches='tight') #%% # #housekeeping ASSIST # #1 days # ndays='1D' # #limit time axis # assistsumplot=assistsum.sel(time=slice(edate-pd.to_timedelta(ndays),edate)) # varsall=[['HBBtopTempassistsum','HBBbottomTempassistsum','HBBapexTempassistsum'],\ # ['ABBtopTempassistsum','ABBbottomTempassistsum','ABBapexTempassistsum'],\ # 'frontEndEnclosureTempassistsum','stepMotorTempassistsum',\ # 'secondPortBlackBodytempassistsum','interferometerTempassistsum',\ # 'ambientAirPowerSupplyTempassistsum','annotatorFtsPlateTempassistsum',\ # 'detectorStirlingCoolerBlockTempassistsum']] # ylabelstr=['$\Delta$ T (°C)','$\Delta$ T (°C)','T (°C)'] # #indicate if difference between all var shall be plotted instead of absolute values 1: yes, 0: no # varsdiff=[1,1,0] # ylims=[[-0.1,0.1],[-0.1,0.1],[]] # fig,axs = plt.subplots(nrows=len(varsall),ncols=1,figsize=(13.5,3*len(varsall)),\ # facecolor='w',edgecolor='w',gridspec_kw={'hspace': 0.07,'wspace': 0.00},\ # tight_layout={'pad':0},sharex=True) # axs=np.atleast_1d(axs) # for i_varsall in range(0,len(varsall)): # vars_=varsall[i_varsall] # if type(axs) is np.ndarray: # ax=axs[i_varsall] # else: # ax=axs # ax.set_prop_cycle(tscycler) # for i_vars in range(0,len(vars_)): # var=vars_[i_vars] # if 'assistsum' in var: # fvar=var[:var.index('assist')] # plotdt=mdates.date2num(assistsumplot['time'].values) # if varsdiff[i_varsall]==1: # varsd_=[v for v in vars_ if not v in var] # for i_varsd in range(len(varsd_)): # vard=varsd_[i_varsd] # fvard=vard[:vard.index('assist')] # plotvar=assistsumplot[fvar]-assistsumplot[fvard] # labelstr=fvar+'-'+fvard # ax.plot(plotdt,plotvar,label=labelstr) # else: # labelstr=fvar # plotvar=assistsumplot[fvar] # ax.plot(plotdt,plotvar,label=labelstr) # ax.grid('on') # try: # ax.set_ylim(ylims[i_varsall]) # except: # print('not ylim defined') # ax.set_ylabel(ylabelstr[i_varsall]) # lines,labels=ax.get_legend_handles_labels() # if len(lines)>8: # ax.legend(fontsize='medium',ncol=int(len(lines)/3.)) # else: # ax.legend(fontsize='medium',ncol=int(len(lines)/2.)) # # if varsdiff[i_varsall]==1: # # varsd_=[v for v in vars_ if not v in var] # # for i_varsd in range(1,len(varsd_)): # # vard=varsd_[i_varsd] # # fvard=vard[:vard.index('assist')] # # plotvar=assistsumplot[fvar]-assistsumplot[fvard] # # ax2=ax.twinx() # # ax2.plot(plotdt,plotvar)