In [1]:
import warnings
warnings.filterwarnings("ignore",category=FutureWarning)
import xarray as xr
import numpy as np
from matplotlib import pyplot as plt
import os
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
warnings.filterwarnings("ignore",message='invalid value encountered in less_equal')
%matplotlib inline
In [2]:
savePath = '/Users/danstechman/GoogleDrive/School/Research/PECAN/Microphysics/plots/vertProfiles_obsOnly'
fType = 'pdf'
noDispSave = True
plotSepSprls = True
plotNt = True
plotTWC = True
plotDmm = True
plotARat = True
plotRE = True
# Define temp bin interval
binIntvl = 1.0
# flights = ['20150617','20150620','20150701','20150702','20150706','20150709']
flights = ['20150620']
In [3]:
tempLim = (22,-18.5)
NtLim = (1e-6,1)
twcLim = (1e-5,2)
DmmLim = (0,2) # 3.5 if all flights
ARatLim = (0.2,0.85)
reLim = (0.01,0.11) # 0.19 if all
In [4]:
# Define temperature bin edges and determine bin midpoints
edgeMin = -19.0 - (binIntvl/2.)
edgeMax = 20.5 + (binIntvl/2.)
edgesTemp = np.arange(edgeMin,edgeMax,binIntvl)
bin_mid = (edgesTemp[0:-1] + edgesTemp[1:])/2
numBins = len(edgesTemp)-1
for flight in flights:
print('Working on {}...'.format(flight))
figSavePath = '{}/{}'.format(savePath,flight)
if not os.path.exists(figSavePath):
os.makedirs(figSavePath)
cipFile = '/Users/danstechman/GoogleDrive/PECAN-Data/mp-data/' + flight + '/' + flight + '_CIPobs-spirals-10s1sAvg.nc'
pecanPrmF = '/Users/danstechman/GoogleDrive/PECAN-Data/' + flight + '_PECANparams.nc'
# Pull out any PECAN parameters
pecanPrms = xr.open_dataset(pecanPrmF,decode_times=False)
mlBotTemp = pecanPrms.mlBotTemp.data
mlTopTemp = pecanPrms.mlTopTemp.data
# Pull out any global variables/attributes from the netcdf file
cipData_root = xr.open_dataset(cipFile)
sprlZone = str(cipData_root.sprlZone.data,'utf-8')
mcsType = str(cipData_root.mcsType.data,'utf-8')
numSprls = len(sprlZone)
# Loop over each spiral for the current flight
for ix in np.arange(0,numSprls):
print('\tAnalyzing Spiral {}'.format(ix+1))
######### Dicitionaries to hold spiral zone data for current spiral
sprlData = {'tempC': [],'rh': [],'Nt': [],'twc': [],'Dmm': [],'ar': [],'re': []}
# Open the group associated with the current spiral
cipData = xr.open_dataset(cipFile,group='spiral_' + str(ix+1))
sprlData['tempC'].append(cipData.tempC_10s.data.tolist())
sprlData['rh'].append(cipData.rh_10s.data.tolist())
sprlData['Nt'].append(cipData.cipNt.data.tolist())
sprlData['twc'].append(cipData.cipTWC_mlr.data.tolist())
sprlData['Dmm'].append(cipData.cipDmm_mlr.data.tolist())
sprlData['ar'].append(cipData.areaRatio_10s.data.tolist())
sprlData['re'].append(cipData.efctvRadius_10s_mlr.data.tolist())
# Pull out all the data for each variable and place within a single list
sprlData['tempC'] = [i for sublist in sprlData['tempC'] for i in sublist]
sprlData['rh'] = [i for sublist in sprlData['rh'] for i in sublist]
sprlData['Nt'] = [i for sublist in sprlData['Nt'] for i in sublist]
sprlData['twc'] = [i for sublist in sprlData['twc'] for i in sublist]
sprlData['Dmm'] = [i for sublist in sprlData['Dmm'] for i in sublist]
sprlData['ar'] = [i for sublist in sprlData['ar'] for i in sublist]
sprlData['re'] = [i for sublist in sprlData['re'] for i in sublist]
# Convert the lists in each dictionary to numpy arrays
sprlData = {key: np.array(val) for key, val in sprlData.items()}
# Change units as desired for any variables
sprlData['re'] = sprlData['re']/1000 # Convert from um to mm
# Set 0's to NaNs in variables where the 0's cause issues
# with the fill plotting (and where a 0 vs. a NaN do not
# make a difference scientifically in my case)
sprlData['Dmm'][sprlData['Dmm'] == 0] = np.nan
sprlData['Nt'][sprlData['Nt'] == 0] = np.nan
sprlData['twc'][sprlData['twc'] == 0] = np.nan
########### Initialize temperature bins and create empty variables for our stats
# Determine which bins each of the temperatures correspond to within each MCS zone
whichBinTemp_sprl = np.digitize(sprlData['tempC'],edgesTemp)
# Define arrays filled with NaNs to hold the min/max/mean/quantiles
# for each variable and at each temperature bin
binRH_mean_sprl,binNt_mean_sprl,binTWC_mean_sprl,binDmm_mean_sprl,binARat_mean_sprl,binRE_mean_sprl, \
= [np.full(numBins,np.nan) for i in range(6)]
########### Bin data by temperature and calculate stats
# Loop through the temperature bins and determine the indices of the
# temperature variable corresponding to temps within said bin
# Then, use these indices to refer to the appropriate values in each of
# our variables of interest
with warnings.catch_warnings():
# Many of our variables have temp bins with all NaNs which will
# throw runtime warnings everytime we try to use nan*math functions
# Here we just tell python to ignore these specific warnings to unclutter
# the output
warnings.filterwarnings('ignore', 'All-NaN (slice|axis) encountered')
warnings.filterwarnings('ignore', 'Mean of empty slice')
for ib in range(0,numBins):
binMatch_sprl = np.squeeze(np.where(whichBinTemp_sprl == ib))
binRH_sprl = sprlData['rh'][binMatch_sprl]
binNt_sprl = sprlData['Nt'][binMatch_sprl]
binTWC_sprl = sprlData['twc'][binMatch_sprl]
binDmm_sprl = sprlData['Dmm'][binMatch_sprl]
binARat_sprl = sprlData['ar'][binMatch_sprl]
binRE_sprl = sprlData['re'][binMatch_sprl]
if np.any(binMatch_sprl):
binRH_mean_sprl[ib] = np.nanmean(binRH_sprl)
binNt_mean_sprl[ib] = np.nanmean(binNt_sprl)
binTWC_mean_sprl[ib] = np.nanmean(binTWC_sprl)
binDmm_mean_sprl[ib] = np.nanmean(binDmm_sprl)
binARat_mean_sprl[ib] = np.nanmean(binARat_sprl)
binRE_mean_sprl[ib] = np.nanmean(binRE_sprl)
########### Single Spiral Plotting ###########
if plotSepSprls:
if plotNt:
fig2, ax2 = plt.subplots(figsize=(16,20))
ax2.axhline(y=mlTopTemp[ix],color='black',linestyle='--',linewidth=3,label='ML Top')
ax2.axhline(y=mlBotTemp[ix],color='red',linestyle='--',linewidth=3,label='ML Bottom')
ax2.plot(sprlData['Nt'],sprlData['tempC'],'bo',markeredgecolor='white',markeredgewidth=1,markersize=12)
ax2.invert_yaxis()
ax2.set_xlim(NtLim)
ax2.set_ylim(tempLim)
ax2.set_xscale('log',nonposx='mask')
ax2.set_xlabel('Total Number Concentration ($cm^{-3}$)',fontsize=24)
ax2.set_ylabel('Temperature ($^{\circ}$C)',fontsize=24)
ax2.tick_params(axis='both', which='major', labelsize=22)
ax2.set_title('{} - $N_t$ - Spiral {}'.format(flight,ix+1),fontsize=26)
ax2.grid(which='both')
ax2.legend(loc='best',fontsize=18)
# Save the output figure
saveStr = '{}/{}_Nt_sprl{:02d}.{}'.format(figSavePath,flight,ix+1,fType)
if noDispSave:
fig2.savefig(saveStr,bbox_inches='tight')
plt.close(fig2)
if plotTWC:
fig3, ax3 = plt.subplots(figsize=(16,20))
ax3.axhline(y=mlTopTemp[ix],color='black',linestyle='--',linewidth=3,label='ML Top')
ax3.axhline(y=mlBotTemp[ix],color='red',linestyle='--',linewidth=3,label='ML Bottom')
ax3.plot(sprlData['twc'],sprlData['tempC'],'bo',markeredgecolor='white',markeredgewidth=1,markersize=12)
ax3.invert_yaxis()
ax3.set_xscale('log',nonposx='mask')
ax3.set_xlim(twcLim)
ax3.set_ylim(tempLim)
ax3.set_xlabel('Total Water Content ($g\ m^{-3}$)',fontsize=24)
ax3.set_ylabel('Temperature ($^{\circ}$C)',fontsize=24)
ax3.tick_params(axis='both', which='major', labelsize=22)
ax3.set_title('{} - TWC (ML Data Excld) - Spiral {}'.format(flight,ix+1),fontsize=26)
ax3.grid(which='both')
ax3.legend(loc='best',fontsize=18)
# Save the output figure
saveStr = '{}/{}_TWC_sprl{:02d}.{}'.format(figSavePath,flight,ix+1,fType)
if noDispSave:
fig3.savefig(saveStr,bbox_inches='tight')
plt.close(fig3)
if plotDmm:
fig4, ax4 = plt.subplots(figsize=(16,20))
ax4.axhline(y=mlTopTemp[ix],color='black',linestyle='--',linewidth=3,label='ML Top')
ax4.axhline(y=mlBotTemp[ix],color='red',linestyle='--',linewidth=3,label='ML Bottom')
ax4.plot(sprlData['Dmm'],sprlData['tempC'],'bo',markeredgecolor='white',markeredgewidth=1,markersize=12)
ax4.invert_yaxis()
ax4.set_xlim(DmmLim)
ax4.set_ylim(tempLim)
ax4.set_xlabel('Median Mass Diameter (mm)',fontsize=24)
ax4.set_ylabel('Temperature ($^{\circ}$C)',fontsize=24)
ax4.tick_params(axis='both', which='major', labelsize=22)
ax4.set_title('{} - $D_{{mm}}$ (ML Data Excld) - Spiral {}'.format(flight,ix+1),fontsize=26)
ax4.xaxis.set_minor_locator(MultipleLocator(0.1))
ax4.grid(which='both')
ax4.legend(loc='best',fontsize=18)
# Save the output figure
saveStr = '{}/{}_Dmm_sprl{:02d}.{}'.format(figSavePath,flight,ix+1,fType)
if noDispSave:
fig4.savefig(saveStr,bbox_inches='tight')
plt.close(fig4)
if plotARat:
fig5, ax5 = plt.subplots(figsize=(16,20))
ax5.axhline(y=mlTopTemp[ix],color='black',linestyle='--',linewidth=3,label='ML Top')
ax5.axhline(y=mlBotTemp[ix],color='red',linestyle='--',linewidth=3,label='ML Bottom')
ax5.plot(sprlData['ar'],sprlData['tempC'],'bo',markeredgecolor='white',markeredgewidth=1,markersize=12)
ax5.invert_yaxis()
ax5.set_xlim(ARatLim)
ax5.set_ylim(tempLim)
ax5.set_xlabel('Area Ratio (%)',fontsize=24)
ax5.set_ylabel('Temperature ($^{\circ}$C)',fontsize=24)
ax5.tick_params(axis='both', which='major', labelsize=22)
ax5.set_title('{} - Area Ratio - Spiral {}'.format(flight,ix+1),fontsize=26)
ax5.xaxis.set_minor_locator(MultipleLocator(0.05))
ax5.grid(which='both')
ax5.legend(loc='best',fontsize=18)
# Save the output figure
saveStr = '{}/{}_ARatio_sprl{:02d}.{}'.format(figSavePath,flight,ix+1,fType)
if noDispSave:
fig5.savefig(saveStr,bbox_inches='tight')
plt.close(fig5)
if plotRE:
fig6, ax6 = plt.subplots(figsize=(16,20))
ax6.axhline(y=mlTopTemp[ix],color='black',linestyle='--',linewidth=3,label='ML Top')
ax6.axhline(y=mlBotTemp[ix],color='red',linestyle='--',linewidth=3,label='ML Bottom')
ax6.plot(sprlData['re'],sprlData['tempC'],'bo',markeredgecolor='white',markeredgewidth=1,markersize=12)
ax6.invert_yaxis()
ax6.set_xlim(reLim)
ax6.set_ylim(tempLim)
ax6.set_xlabel('Effective Radius (mm)',fontsize=24)
ax6.set_ylabel('Temperature ($^{\circ}$C)',fontsize=24)
ax6.tick_params(axis='both', which='major', labelsize=22)
ax6.set_title('{} - $R_e$ (ML Data Excld) - Spiral {}'.format(flight,ix+1),fontsize=26)
ax6.xaxis.set_minor_locator(MultipleLocator(0.01))
ax6.grid(which='both')
ax6.legend(loc='best',fontsize=18)
# Save the output figure
saveStr = '{}/{}_RE_sprl{:02d}.{}'.format(figSavePath,flight,ix+1,fType)
if noDispSave:
fig6.savefig(saveStr,bbox_inches='tight')
plt.close(fig6)