# Figure 5 of paper "A locally smoothed terrain-following coordinate" # submitted to JAMES # Author: Stephanie Westerhuis # Date: February 1, 2021 import matplotlib import numpy as np import xarray as xr import matplotlib as mpl import matplotlib.pyplot as plt import pandas as pd import datetime as dt import sys from scipy import interpolate, stats import math from pathlib import Path import cartopy import matplotlib.image as mpimg from matplotlib.patches import Patch from matplotlib.lines import Line2D import seaborn as sns from scipy import spatial from variables import temp, windspeed, qv, qc # compare vertical profiles of # temparature, wind, specific humidity, liquid water content def get_latlon_kdtree(lats, lons): return spatial.KDTree(np.array(list(zip( lons.ravel(), lats.ravel())))) def get_ind_closest_to_latlon(lat, lon, lats, lons): kdtree = get_latlon_kdtree(lats, lons) x_shp = lons.shape[0] y_shp = lons.shape[1] data_index = np.arange(0, x_shp*y_shp).reshape(x_shp, y_shp) distance, i = kdtree.query((lon, lat)) ind_lat_point, ind_lon_point = np.argwhere(data_index == i)[0] return ind_lat_point, ind_lon_point def calc_alt(hhl): # altitude: height of full levels return hhl[1:,:,:]+(hhl[:-1,:,:]-hhl[1:,:,:])/2 def calc_qv_from_td(td, p): """ td: dewpoint temp in °C p: pressure in hPa """ e = 6.112*np.exp((17.67*td)/(td + 243.5)) qv = (0.622 * e)/(p - (0.378 * e)) * 1000 # to get g/kg # equivalent formula would be: # from calculating spec. hum. from dew point temperatures github-discussion # after eq. 4.24 in Practical Meteorology from Stull # td=td+273.15 # qv = (622*6.113*np.exp(5423*(td-273.15)/(td*273.15)))/p return qv # files lts = ['0012','0015','0018','0021','0100'] d1 = 'data/fcst_sleve' # SLEVE d6 = 'data/fcst_losvec' # LOSVEC fs1 = [Path(d1, f'lfff{lt}0000.nc') for lt in lts] fs6 = [Path(d6, f'lfff{lt}0000.nc') for lt in lts] ds1 = xr.open_mfdataset(fs1, combine='by_coords') ds6 = xr.open_mfdataset(fs6, combine='by_coords') # location payerne (lat_pay, lon_pay) = (46.81303, 6.94302) # sleve c1 = xr.open_dataset(Path('data','const_sleve.nc')).isel(time=0) lats1 = c1.lat_1.values lons1 = c1.lon_1.values lat1, lon1 = get_ind_closest_to_latlon(lat_pay,lon_pay,lats1,lons1) hhl1 = c1.HHL.values alt1 = calc_alt(hhl1)[:,lat1,lon1] # losvec c6 = xr.open_dataset(Path('data','const_losvec.nc')).isel(time=0) lats6 = c6.lat.values lons6 = c6.lon.values lat6, lon6 = get_ind_closest_to_latlon(lat_pay,lon_pay,lats6,lons6) hhl6 = c6.HHL.values alt6 = calc_alt(hhl6)[:,lat6,lon6] # colors c1 = 'yellowgreen' c6 = 'mediumblue' nrows = 4 ncols = 3 fig, axs = plt.subplots(nrows=nrows, ncols=ncols, figsize=(7.25,10), sharey=True, dpi=300) fig.subplots_adjust(left=0.1, bottom=0.03, right=0.97, top=0.97, wspace=0.15, hspace=0.3) for (row, var) in [(0,temp),(1,windspeed),(2,qv)]: if var.shortname == 'windspeed': uwind1 = ds1.U.values[:,:,lat1,lon1] vwind1 = ds1.V.values[:,:,lat1,lon1] v1 = np.sqrt(uwind1**2 + vwind1**2) uwind6 = ds6.U.values[:,:,lat6,lon6] vwind6 = ds6.V.values[:,:,lat6,lon6] v6 = np.sqrt(uwind6**2 + vwind6**2) else: v1 = ds1[var.cosmoname].values[:,:,lat1,lon1] \ * var.mult_factor v6 = ds6[var.cosmoname].values[:,:,lat6,lon6] \ * var.mult_factor # 00 UTC axs[row,0].plot(v1[0,:,], alt1, color = c1, linestyle = '-', label = 'SLEVE', ) axs[row,0].plot(v6[0,:,], alt6, color = c6, linestyle = '--', label = 'LOSVEC', ) # 06 UTC axs[row,1].plot(v1[2,:,], alt1, color = c1, linestyle = '-', label = 'SLEVE', ) axs[row,1].plot(v6[2,:,], alt6, color = c6, linestyle = '--', label = 'LOSVEC', ) # 12 UTC axs[row,2].plot(v1[4,:,], alt1, color = c1, linestyle = '-', label = 'SLEVE', ) axs[row,2].plot(v6[4,:,], alt6, color = c6, linestyle = '--', label = 'LOSVEC', ) # qc plots var = qc row = 3 v1 = ds1[var.cosmoname].values[:,:,lat1,lon1] \ * var.mult_factor v6 = ds6[var.cosmoname].values[:,:,lat6,lon6] \ * var.mult_factor # 03 UTC axs[row,0].plot(v1[1,:,], alt1, color = c1, linestyle = '-', label = 'SLEVE', ) axs[row,0].plot(v6[1,:,], alt6, color = c6, linestyle = '--', label = 'LOSVEC', ) # 06 UTC axs[row,1].plot(v1[2,:,], alt1, color = c1, linestyle = '-', label = 'SLEVE', ) axs[row,1].plot(v6[2,:,], alt6, color = c6, linestyle = '--', label = 'LOSVEC', ) # 09 UTC axs[row,2].plot(v1[3,:,], alt1, color = c1, linestyle = '-', label = 'SLEVE', ) axs[row,2].plot(v6[3,:,], alt6, color = c6, linestyle = '--', label = 'LOSVEC', ) ## add radiosounding rs00 = 'data/2016122900_pay_sounding.txt' rs12 = 'data/2016122912_pay_sounding.txt' d00 = pd.read_csv(rs00, sep=";",na_values=['na']) d12 = pd.read_csv(rs12, sep=";",na_values=['na']) # temp for (col, rs) in [(0,d00),(2,d12)]: tt = rs.zrettts0+273 hh = rs.zreppzs0 axs[0,col].plot(tt,hh, color = 'k', linestyle = ':', label = 'RS', ) # windspeed for (col, rs) in [(0,d00),(2,d12)]: ww = rs.zklfffs0 hh = rs.zreppzs0 axs[1,col].plot(ww,hh, color = 'k', linestyle = ':', label = 'RS', ) # specific humidity for (col, rs) in [(0,d00),(2,d12)]: td = rs.zretdes0 # in °C p = rs.zreppps0 # in hPa sh = calc_qv_from_td(td, p) hh = rs.zreppzs0 axs[2,col].plot(sh,hh, color = 'k', linestyle = ':', label = 'RS', ) # legend fig.axes[2].legend() for ax in fig.axes: ax.set_ylim(490,1600) # temp for ax in fig.axes[0:3]: ax.set_xlim(267,282) # wind for ax in fig.axes[3:6]: ax.set_xlim(0,15) # qv for ax in fig.axes[6:9]: ax.set_xlim(0.4,3.8) # qc for ax in fig.axes[9:12]: ax.set_xlim(-0.01,0.28) axs[0,0].set_title('Temp. [K]: 00 UTC', loc='left') axs[0,0].set_ylabel('Elevation [m asl]') axs[1,0].set_title('Wind [m/s]: 00 UTC', loc='left') axs[1,0].set_ylabel('Elevation [m asl]') axs[2,0].set_title('Hum. [g/kg]: 00 UTC', loc='left') axs[2,0].set_ylabel('Elevation [m asl]') axs[3,0].set_title('LWC [g/kg]: 03 UTC', loc='left') axs[3,0].set_ylabel('Elevation [m asl]') axs[0,1].set_title('06 UTC') axs[0,2].set_title('12 UTC') axs[3,2].set_title('09 UTC') # label figures labels = ['a','b','c','d','e','f','g','h', 'i','j','k','l'] for ii, ax in enumerate(fig.axes): ax.text(0.94,0.025,labels[ii], horizontalalignment='left', transform=ax.transAxes) # save out_name = '5_profiles' fig.savefig(f'plots/{out_name}.png')