# -*- coding: utf-8 -*- """ Created on Wed Nov 3 16:10:20 2021 @author: Sreenath """ #%% clear environment try: from IPython import get_ipython get_ipython().magic('clear') get_ipython().magic('reset -f') except: pass #%% libraries import numpy as np import xarray as xr import matplotlib.pyplot as plt import pandas as pd import shutil, os import matplotlib as mpl import seaborn as sns #%% plot parameters mpl.rcParams['axes.linewidth'] = 2.0 #set the value globally plt.rc('font', family='serif',size = 16) # controls default text layout plt.rc('axes', titlesize=16) # fontsize of the axes title plt.rc('axes', labelsize=16) # fontsize of the x and y labels plt.rc('xtick', labelsize=10) # fontsize of the tick labels plt.rc('ytick', labelsize=16) # fontsize of the tick labels plt.rc('legend', fontsize=16) # legend fontsize plt.rc('figure', titlesize=16) # fontsize of the figure title #plt.minorticks_on() plt.rc('text', usetex=False) sns.set(font = "Cambria",font_scale=1.5,style="whitegrid") #sns.set(font = "Cambria",font_scale=1.5,style="white") #mpl.rc_file_defaults() #undos seaborn effects #custom fontsize font = 12 #%% data folder source = r'C:\Users\Sreenath\Documents\palm\Cheyenne\realistic_IOP03\OUTPUT\run_1' os.chdir(source) #get list of all the folders folder_list = os.listdir(source) print(folder_list) #data file name to concatenate input_fname = 'DATA_1D_PR_NETCDF' domain = '_N02' pr_data = input_fname + domain os.chdir(source) #combine profile data pr_list = [] for i, folder in enumerate(folder_list): #open data file inside the folder. ds = xr.open_dataset(source+ '\\' + folder + '\\' + pr_data) #append the data file name to the list pr_list.append(ds) #concatenate all data files in the list pr_ds = xr.concat(pr_list, dim='time') print('done!') #%% zu_simulated = np.asarray(pr_ds['zu']) zs_simulated = np.asarray(pr_ds['zs']) zv_simulated = np.asarray(pr_ds['zu']) ztheta_simulated = np.asarray(pr_ds['ztheta']) zq_simulated = np.asarray(pr_ds['zq']) zwtheta_simulated = np.asarray(pr_ds['zwtheta']) #%% plot and save simulated profiles os.chdir(r'C:\Users\Sreenath\Documents\palm\Cheyenne\realistic_IOP03\plots\run.1') time_range = range(1,96,12) #time_range = range(96,193,12) day = '_day_1' fig = plt.figure(figsize=(8,8)) ax = fig.subplots(1,1) #ax = plt.subplot(111) #for i in range(1,len(pr_ds.time),6): #for i in range(96,len(pr_ds.time),12): #second day hours for i in time_range: #first day hours ax.plot(np.asarray(pr_ds.theta.sel(time=pr_ds.time[i])),ztheta_simulated,'-',label=str(pr_ds.time[i].data/(3600*1E9)))#1E9*3600 to convert to hours #ax.set_ylim([0., 100.]) ax.set_ylabel('z (m)',fontsize=18) ax.set_xlabel('theta simulated (K)',fontsize=18) fig.legend(loc=7,fontsize=12) fig.tight_layout() fig.subplots_adjust(right=0.75) fname = 'theta' + domain + day plt.savefig(fname,dpi=300) fig = plt.figure(figsize=(8,8)) ax = fig.subplots(1,1) for i in time_range: ax.plot(np.asarray(pr_ds.w.sel(time=pr_ds.time[i])),zwtheta_simulated,'-',label=str(pr_ds.time[i].data/(3600*1E9)))#1E9*3600 to convert to hours, #ax.set_ylim([0., 100.]) ax.set_xlabel('w (m/s)',fontsize=18) ax.set_ylabel('z (m)',fontsize=18) fig.legend(loc=7,fontsize=12) fig.tight_layout() fig.subplots_adjust(right=0.75) fname = 'w' + domain + day plt.savefig(fname,dpi=300) fig = plt.figure(figsize=(8,8)) ax = fig.subplots(1,1) for i in time_range: ax.plot(np.asarray(pr_ds.wtheta.sel(time=pr_ds.time[i])),zwtheta_simulated,'-',label=str(pr_ds.time[i].data/(3600*1E9)))#1E9*3600 to convert to hours ax.set_ylim([0., 2200.]) ax.set_xlim([-100., 300.]) ax.set_xlabel('wtheta (W/m2)',fontsize=18) ax.set_ylabel('z (m)',fontsize=18) fig.legend(loc=7,fontsize=12) fig.tight_layout() fig.subplots_adjust(right=0.75) fname = 'wtheta' + domain + day plt.savefig(fname,dpi=300) # fig = plt.figure(figsize=(8,8)) # ax = fig.subplots(1,1) # for i in range(1,len(pr_ds.time),12): # wind = (pr_ds.u.sel(time=pr_ds.time[i])**2+pr_ds.v.sel(time=pr_ds.time[i])**2)**(0.5) # ax.plot(np.asarray(wind[0]),zu_simulated,'-',label=str(pr_ds.time[i].data/(1E9*3600))) # ax.set_ylim([0., 1500.]) # ax.set_xlabel('horizontal wind (m/s)',fontsize=18) # ax.set_ylabel('z (m)',fontsize=18) # fig.legend(fontsize=12) fig = plt.figure(figsize=(8,8)) (ax1,ax2) = fig.subplots(1,2,sharex=False, sharey=True) #for i in range(0,len(pr_ds.time),6): for i in time_range: ax1.plot(np.asarray(pr_ds.u.sel(time=pr_ds.time[i])),zu_simulated,'-') #ax1.set_ylim([0., 3000.]) ax2.plot(np.asarray(pr_ds.v.sel(time=pr_ds.time[i])),zu_simulated,'-',label=str(pr_ds.time[i].data/(3600*1E9)))#1E9*3600) #ax2.set_ylim([0., 3000.]) #ax1.set_xlim([-2., 15.]) #ax1.set_ylim([0., 1500.]) #ax2.set_ylim([0., 1500.]) ax1.set_xlabel('u (m/s)',fontsize=18) ax2.set_xlabel('v (m/s)',fontsize=18) ax1.set_ylabel('z (m)',fontsize=18) fig.legend(loc=7,fontsize=12) fig.tight_layout() fig.subplots_adjust(right=0.75) fname = 'u_v' + domain + day plt.savefig(fname,dpi=300) fig = plt.figure(figsize=(8,8)) ax = fig.subplots(1,1) for i in time_range: ax.plot(np.asarray(pr_ds.q.sel(time=pr_ds.time[i])),ztheta_simulated,'-',label=str(pr_ds.time[i].data/(3600*1E9))) #ax.set_ylim([0., 100.]) ax.set_ylabel('z (m)',fontsize=18) ax.set_xlabel('wv mixing ratio (kg/kg)',fontsize=18) fig.legend(loc=7,fontsize=12) fig.tight_layout() fig.subplots_adjust(right=0.75) fname = 'q' + domain + day plt.savefig(fname,dpi=300) fig = plt.figure(figsize=(8,8)) ax = fig.subplots(1,1) for i in time_range: ax.plot(np.asarray(pr_ds.wq.sel(time=pr_ds.time[i])),zwtheta_simulated,'-',label=str(pr_ds.time[i].data/(3600*1E9)))#1E9*3600 to convert to hours ax.set_ylim([0., 2200.]) ax.set_xlabel('wq (W/m2)',fontsize=18) ax.set_ylabel('z (m)',fontsize=18) fig.legend(loc=7,fontsize=12) fig.tight_layout() fig.subplots_adjust(right=0.75) fname = 'wq' + domain + day plt.savefig(fname,dpi=300) fig = plt.figure(figsize=(8,8)) ax = fig.subplots(1,1) #for i in range(1,len(pr_ds.time),6): for i in time_range: #first 4 hours ax.plot(np.asarray(pr_ds.s.sel(time=pr_ds.time[i])),zs_simulated,'-',label=str(pr_ds.time[i].data/(3600*1E9)))#1E9*3600 to convert to hours ax.set_ylim([0., 1500.]) ax.set_ylabel('z (m)',fontsize=18) ax.set_xlabel('CO2 simulated (mg/Kg)',fontsize=18) fig.legend(loc=7,fontsize=12) fig.tight_layout() fig.subplots_adjust(right=0.75) fname = 'co2' + domain + day plt.savefig(fname,dpi=300) #%%#plotting IOP sonde profiles source = r'C:\Users\Sreenath\Documents\palm\Cheyenne\ches_forcings\ISS_sonde' folder = '20190924' os.chdir(source + '\\' + folder ) #pr_sonde_data = 'CHEESEHEAD_ISS1_20190924_175905_RS41_isf.nc' #get list of all the files file_list = os.listdir(source + '\\' + folder) timelist = ['05:47','09:15','12:59','16:45'] fig = plt.figure(figsize=(8,8)) ax = fig.subplots(1,1) for i, file in enumerate(file_list): pr_sonde_data = file pr_sonde_ds = xr.open_dataset(file) #ax.plot(np.asarray(pr_sonde_ds['mr']),np.asarray(pr_sonde_ds['alt']),'-',) #ax.plot(np.asarray(pr_sonde_ds['theta']+273.15),np.asarray(pr_sonde_ds['alt']-462),'-',label=str(timelist[i])) ax.plot(np.asarray(pr_sonde_ds['mr']),np.asarray(pr_sonde_ds['alt']-462),'-',label=(timelist[i])) ax.set_ylim([0, 1750.]) #ax.set_xlim([280, 310]) ax.set_ylabel('z (m)',fontsize=18) ax.set_xlabel('mr sonde (g/Kg)',fontsize=18) fig.legend(fontsize=12) #%% compare simulated N02 profiles with the IOP sonde profiles time_range = range(1,96,12) #time_range = range(96,193,12) day = '_day_1' simulated_time_index = [23,37,52,67] simulated_time = ['05:45','09:15','13:00','16:45'] fig = plt.figure(figsize=(8,8)) ax = fig.subplots(1,1) #for i in range(1,len(pr_ds.time),6): #for i in range(96,len(pr_ds.time),12): #second day hours for i in range(0,4): #first day hours # ax.plot(np.asarray(pr_ds.theta.sel(time=pr_ds.time[i])),ztheta_simulated,'-',label=str(pr_ds.time[i].data/(3600*1E9)))#1E9*3600 to convert to hours ax.plot(np.asarray(pr_ds.theta.sel(time=pr_ds.time[simulated_time_index[i]])),ztheta_simulated,'--',label=simulated_time[i])#1E9*3600 to convert to hours ax.set_ylim([0, 1750.]) #ax.set_xlim([290, 310]) ax.set_ylabel('z (m)',fontsize=18) ax.set_xlabel('theta simulated (K)',fontsize=18) fig.legend(loc=7,fontsize=12) fig.tight_layout() fig.subplots_adjust(right=0.75) #fname = 'theta' + domain + day #plt.savefig(fname,dpi=300) #%% input profile data and time, color lists source = r'C:\Users\Sreenath\Documents\palm\Cheyenne\ches_forcings\ISS_sonde' #day_1 # folder = '20190924' # timelist = ['05:47','09:15','12:59','16:45'] #day_2 folder = '20190925' timelist = ['05:50','09:13','12:57','16:45'] os.chdir(source + '\\' + folder ) #pr_sonde_data = 'CHEESEHEAD_ISS1_20190924_175905_RS41_isf.nc' #get list of all the files file_list = os.listdir(source + '\\' + folder) #day_1 #simulated_time_index = [23,37,52,67] simulated_time_index = [119,133,148,163] simulated_time = ['05:45','09:15','13:00','16:45'] color_list = ['blue','orange','green','red'] #%%#theta, q profiles fig = plt.figure(figsize=(8,8)) ax = fig.subplots(1,1) for i, file in enumerate(file_list): pr_sonde_data = file pr_sonde_ds = xr.open_dataset(file) ax.plot(np.asarray(pr_sonde_ds['theta']+273.15),np.asarray(pr_sonde_ds['alt']-462),'-',label=str(timelist[i]),color=color_list[i]) #ax.plot(np.asarray(pr_sonde_ds['mr']),np.asarray(pr_sonde_ds['alt']-462),'-',label=str(timelist[i]),color=color_list[i]) for i in range(0,4): #first day hours ax.plot(np.asarray(pr_ds.theta.sel(time=pr_ds.time[simulated_time_index[i]])),ztheta_simulated,'--',label=simulated_time[i],color = color_list[i]) #ax.plot(np.asarray(pr_ds.q.sel(time=pr_ds.time[simulated_time_index[i]])*1000),ztheta_simulated,'--',label=simulated_time[i],color = color_list[i]) #ax.set_xlabel('Theta (K)',fontsize=18) ax.set_ylim([0, 1750.]) #ax.set_xlim([280, 310]) ax.set_xlim([284, 301]) ax.set_ylabel('z (m)',fontsize=18) #ax.set_xlabel('Mixing ratio (g/Kg)',fontsize=18) ax.set_xlabel('Theta (K)',fontsize=18) fig.legend(fontsize=12) #%%# u v wind profiles fig = plt.figure(figsize=(8,8)) (ax1,ax2) = fig.subplots(1,2,sharex=False, sharey=True) #sonde profiles for i, file in enumerate(file_list): pr_sonde_data = file pr_sonde_ds = xr.open_dataset(file) ax1.plot(np.asarray(pr_sonde_ds['u_wind']),np.asarray(pr_sonde_ds['alt']-462),'-',label=str(timelist[i]),color=color_list[i]) #ax1.set_ylim([0., 3000.]) ax2.plot(np.asarray(pr_sonde_ds['v_wind']),np.asarray(pr_sonde_ds['alt']-462),'-',label=str(timelist[i]),color=color_list[i]) #ax2.set_ylim([0., 3000.]) for i in range(0,4): #first day hours ax1.plot(np.asarray(pr_ds.u.sel(time=pr_ds.time[simulated_time_index[i]])),zu_simulated,'--',label=simulated_time[i],color = color_list[i]) #ax.set_xlabel('Theta (K)',fontsize=18) ax2.plot(np.asarray(pr_ds.v.sel(time=pr_ds.time[simulated_time_index[i]])),zv_simulated,'--',label=simulated_time[i],color = color_list[i]) #ax1.set_xlim([-2., 15.]) ax1.set_ylim([0., 1750.]) ax1.set_xlim([-5., 20.]) ax2.set_ylim([0., 1750.]) ax2.set_xlim([-10., 5.]) ax1.set_xlabel('u (m/s)',fontsize=18) ax2.set_xlabel('v (m/s)',fontsize=18) ax1.set_ylabel('z (m)',fontsize=18) ax2.legend(fontsize=12) #%%w wind #%%# u v wind profiles fig = plt.figure(figsize=(8,8)) ax = fig.subplots(1,1) #sonde profiles for i, file in enumerate(file_list): pr_sonde_data = file pr_sonde_ds = xr.open_dataset(file) ax.plot(np.asarray(pr_sonde_ds['dz']),np.asarray(pr_sonde_ds['alt']-462),'-',label=str(timelist[i]),color=color_list[i]) for i in range(0,4): #first day hours ax.plot(np.asarray(pr_ds.w.sel(time=pr_ds.time[simulated_time_index[i]])),zwtheta_simulated,'--',label=simulated_time[i],color = color_list[i]) #ax1.set_xlim([-2., 15.]) ax.set_ylim([0., 1750.]) #ax.set_xlim([-5., 16.]) ax.set_xlabel('w (m/s)',fontsize=18) ax.set_ylabel('z (m)',fontsize=18) fig.legend(fontsize=12)