<jupyter_start><jupyter_text>27.02.2018 -- Email from Steve: also be interesting:    -  calculate the total snow water path (SWP, integrated SWC throughout the depth of the cloud)$\rightarrow$ could then compare the model SWP with the retrieval SWP to see if we are in the right ball parki.e. see if the models are producing the right amount of snow to the first order.  This also may give us further insight into how to deal with the units.<jupyter_code>import sys
sys.path.append('/Volumes/SANDISK128/Documents/Thesis/Python/')
sys.path.append('/Volumes/SANDISK128/Documents/Thesis/Python/weather_mast/')
sys.path.append('/Volumes/SANDISK128/Documents/Thesis/Python/Retrieval/')
import netCDF4
import fill_values as fv
import datetime
import calc_date as cd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import save_fig as sF
import createFolder as cF
from matplotlib.gridspec import GridSpec
import plt_ce_image as im
import read_MRR as pMRR

from scipy.integrate import simps
year = '2016'
month = '12'
#t = ['26']
#t = ['17', '18', '19', '20', '21', '22', '23', '24', '25', '26', '27']
t = ['20', '21', '22', '23', '24', '25', '26', '27']
tid = '00'
station = 'Haukeliseter'


### Figures 
savefig = 0
figdir = '../../Figures/Retrieval_MEPS/snow_graupel_cloudice_SWP_THICKNESS'
cF.createFolder('%s/' %(figdir))
form = 'png'
################################################################
### MEPS
level = 'sfc'
nc_dir_sfc = '../../Data/MEPS/%s/%s_%s' %(station,level,tid)
level = 'ml'
nc_dir_ml = '../../Data/MEPS/%s/%s_%s' %(station,level,tid)
# snow fall amount
var_name_snow = 'snowfall_amount_ml'
snow_dir = '%s/%s' %(nc_dir_ml,var_name_snow)
# graupel fall amount
var_name_graupel = 'graupelfall_amount_ml'
graupel_dir = '%s/%s' %(nc_dir_ml,var_name_graupel)
# atmosphere_cloud_ice_content_ml
var_name_cl_ice = 'atmosphere_cloud_ice_content_ml'
cl_ice_dir = '%s/%s' %(nc_dir_ml,var_name_cl_ice)

################################################################
### MRR
MRR_dir = '../../Data/MRR/processed_MRR'

################################################################
### Retrieval 
nc_dir_retrieval = '../../Data/Retrieved_SWC'

################################################################
snow_file         = dict()
graupel_file      = dict()
cl_ice_file       = dict()

fn_snow           = dict()
fn_graupel        = dict()
fn_cl_ice         = dict()

time_ml           = dict()
snow_amount_ml    = dict()
graupel_amount_ml = dict()
cl_ice_ml         = dict()

thickness_ml      = dict()
ice_amount_ml     = dict()
ice_amount_ml_TH  = dict()
time              = dict()
height_ml         = dict()
height            = dict()
SWP_model         = dict()
for day in t:
## MRR ######  
    fnMRR = netCDF4.Dataset('%s/VMRR_%s%s%s.nc' %(MRR_dir,year,month,day) ,'r')
    time_MRR   = fnMRR.variables['time'][:]
    height_MRR = fnMRR.variables['height'][:]

    Ze = pMRR.read_and_mask(fnMRR,'Ze', np.nan)         # vertical Ze profile for retrieval

## RETRIEVAL ######
    Retrieval_file = '%s/out_swc_b6_%s%s%s.nc' %(nc_dir_retrieval, year, month, day)
    fn = netCDF4.Dataset(Retrieval_file)

    n_profiles = fn.variables['n_profiles'][:]
    n_layers = fn.variables['n_layers'][:]
    retrieved_snowfall_amount = fn.variables['retrieved_snowfall_amount'][:]
    SWC = fn.variables['SWC'][:]
## MEPS ######
    for ens_memb in range(0,10):
        snow_file[ens_memb]    = '%s/%s%s%s_%s_%s.nc' %(snow_dir,year,month,day,tid,ens_memb)
        graupel_file[ens_memb] = '%s/%s%s%s_%s_%s.nc' %(graupel_dir,year,month,day,tid,ens_memb)
        cl_ice_file[ens_memb]  = '%s/%s%s%s_%s_%s.nc' %(cl_ice_dir,year,month,day,tid,ens_memb)
    
        fn_snow[ens_memb]      = netCDF4.Dataset(snow_file[ens_memb])
        fn_graupel[ens_memb]   = netCDF4.Dataset(graupel_file[ens_memb]) 
        fn_cl_ice[ens_memb]    = netCDF4.Dataset(cl_ice_file[ens_memb])
# Read in time to know initialisation time (plotting)    
        time_ml[ens_memb] = fn_snow[ens_memb].variables['time']
        time_ml[ens_memb] = fv.fill_nan(time_ml[ens_memb][:])
        

        ini_day = datetime.datetime.utcfromtimestamp(time_ml[0][0]).day   # day of initialisation
        hh = datetime.datetime.utcfromtimestamp(time_ml[0][0]).hour       # first hour of initialisation?
        calday, calmon = cd.get_dayname(year, month, ini_day)
# Read in the variable name which should be plotted  
        snow_amount_ml[ens_memb]    = fn_snow[ens_memb].variables[var_name_snow]
        graupel_amount_ml[ens_memb] = fn_graupel[ens_memb].variables[var_name_graupel]
        cl_ice_ml[ens_memb]         = fn_cl_ice[ens_memb].variables[var_name_cl_ice]


# create an array with time of the shape of the variable
        lead_time   = np.arange(0,snow_amount_ml[ens_memb].shape[0])
        model_level = np.arange(0,snow_amount_ml[ens_memb].shape[1])
        td,Y        = np.meshgrid(lead_time,model_level)

# substitute missing values with nan       
        snow_amount_ml[ens_memb]    = fv.fill_nan(snow_amount_ml[ens_memb][:])
        graupel_amount_ml[ens_memb] = fv.fill_nan(graupel_amount_ml[ens_memb][:])
        cl_ice_ml[ens_memb]         = fv.fill_nan(cl_ice_ml[ens_memb][:])
# get ice amount (snow+graupel+cloud ice)
        ice_amount_ml[ens_memb]     = snow_amount_ml[ens_memb] + graupel_amount_ml[ens_memb] + cl_ice_ml[ens_memb]
    
# layer thickness to calculate height
        thickness_ml[ens_memb]      = fn_snow[ens_memb].variables['layer_thickness']
        thickness_ml[ens_memb]      = fv.fill_nan(thickness_ml[ens_memb][:])
    
        thickness = thickness_ml[ens_memb]
        h_above = []
        h_above.append(np.zeros(thickness[:,:].shape[0]))
        th_arr = np.nansum([np.asarray(h_above)[:,0], thickness[:,0]], axis = 0)
        h_above.append(th_arr)
        for i in range(2,thickness[:,:].shape[1]):
            th_arr = np.nansum([h_above[i-1], thickness[:,i-1]], axis = 0)
            h_above.append(th_arr)
        height_ml[ens_memb] = np.transpose(h_above)
        
        
# get only valuable values not nan    
        time_ml[ens_memb]       = (time_ml[ens_memb][~np.isnan(time_ml[ens_memb])])
        time[ens_memb]          = (td[:,~np.isnan(ice_amount_ml[ens_memb]).any(axis=1)])
        height[ens_memb]        = (height_ml[ens_memb][~np.isnan(ice_amount_ml[ens_memb]).any(axis=1),:])
        thickness_ml[ens_memb]  = (thickness_ml[ens_memb][~np.isnan(ice_amount_ml[ens_memb]).any(axis=1),:])  
        ice_amount_ml[ens_memb] = (ice_amount_ml[ens_memb][~np.isnan(ice_amount_ml[ens_memb]).any(axis=1)])

# convert the snow amount from [kg/m^2] into [g/m^2] by multiply with 1000
        ice_amount_ml[ens_memb] = (ice_amount_ml[ens_memb]*1000.)
        

# convert the snow amount from [g/m^2] into [g/m^3] by multiplying with the thickness
        
        ice_amount_ml_TH[ens_memb] = (1/thickness_ml[ens_memb]) * ice_amount_ml[ens_memb]
    

        SWP_ml = []
        for i in range(0,ice_amount_ml_TH[ens_memb][:,:].shape[0]):
            SWP_ml.append(simps(y = ice_amount_ml_TH[ens_memb][i,np.where(height[ens_memb][i,:] <= 3000)], 
                            x = height[ens_memb][i,np.where(height[ens_memb][i,:] < 3000)]))
#            SWP_ml.append(simps(y = ice_amount_ml[ens_memb][i,np.where(height[ens_memb][i,:] <= 3000)], 
 #                           x = height[ens_memb][i,np.where(height[ens_memb][i,:] < 3000)]))
        SWP_model[ens_memb] = SWP_ml

        h1, h2 = np.where(height[ens_memb][:,:] > 3000)
        ice_amount_ml[ens_memb][h1[:],h2[:]] = np.nan
        ice_amount_ml_TH[ens_memb][h1[:],h2[:]] = np.nan
        time[ens_memb][h2[:],h1[:]]          = -9999
        height[ens_memb][h1[:],h2[:]]        = np.nan
    

#    for ens_memb in range(0,10):
 #       if len(ice_amount_ml_TH[ens_memb]) == 0:
  #          continue
   #     else:
    #        print('%s%s%s_%s' %(year,month,day,tid))
     #       print('ice_amount_ml_TH',ens_memb,np.nanmax(ice_amount_ml_TH[ens_memb]))#,
      #     #   'ice_cont_ml', ens_memb, np.nanmax(ice_cont_ml[ens_memb]))

    SWP_ret = []
    SWP_ret = simps(y = SWC[:,:], x = height_MRR, axis = 0)

#    print('SWP_ret',np.nanmax(SWP_ret))

#    for ens_memb in range(0,10):
 #       if len(SWP_model[ens_memb]) == 0:
  #          continue
   #     else:
    #        print('SWP_model',ens_memb,np.nanmax(SWP_model[ens_memb]))
###################################################################################
# for plotting, only values larger than 0.01
    a,b = np.where(retrieved_snowfall_amount[:,:] < 0.01)
    SWC[a[:],b[:]] = np.nan
    retrieved_snowfall_amount[a[:],b[:]] = np.nan
    np.transpose(Ze)[a[:],b[:]] = np.nan


#a2,b2 = np.where(snow_amount_ml[0][:,:] < 0.001)
#snow_amount_ml[0][a2[:],b2[:]] = np.nan
#    for ens_memb in range(0,10):
 #       if len(ice_amount_ml[ens_memb]) == 0:
  #          continue
   #     else:
    #        a2,b2 = np.where(ice_amount_ml[ens_memb][:,:] < 0.01)    ## 0.001? retrieved_snowfall_amount = iwcpsd*0.85*e-3*3600*24
     #       ice_amount_ml[ens_memb][a2[:],b2[:]] = np.nan
#     #       ice_amount_ml_TH[ens_memb][a2[:],b2[:]] = np.nan

########################################
### PLOTTING ###########
    fig_name = '%s%s%s.%s' %(year,month,ini_day,form)
    fig = plt.figure(figsize=(18.,30))
# title 

#    fig.suptitle('%s, %s %s %s' %(calday, day, calmon, year), y=1.01,
 #           color =im.date_blue, fontsize=24)
    gs = GridSpec(5,2)
## 1st subplot --- REFLECTIVITY
    ax0 = plt.subplot(gs[0,:])
    
    levels = np.arange(-10,30.1,0.1)
#    levels = np.arange(0,75.1,0.1)
#    if day == '26':
 #       im.plt_ce_image(fig,ax0,time_MRR[:-1],height_MRR,Ze,levels,
  #          v_min=0.0, v_max=75., xmax=time_MRR.max(), ymax=height_MRR.max(), 
   #         xticks=im.xticks1, yticks=im.yticks1,
    #        cb_ticks=levels[::150],var_label='P [mm day$^{-1}$]')
#    else:
    im.plt_ce_image(fig,ax0,time_MRR, height_MRR, np.transpose(Ze),levels,
            v_min=0.0, v_max=30., xmax=time_MRR.max(), ymax=height_MRR.max(), 
            xticks=im.xticks1, yticks=im.yticks1,
            cb_ticks=levels[::50],var_label='Ze [dBZ]')
    ax0.text(0.02,0.95, '%s, %s %s %s' %(calday, day, calmon, year),
            verticalalignment ='top', horizontalalignment='left',
             transform = ax0.transAxes,
             color = im.date_blue, fontsize = im.fontsize,
             bbox={'facecolor':'white','alpha':1.,'pad':10})
    
################
## 2nd subplot --- SNOW WATER CONTENT RETRIEVAL
    ax1 = plt.subplot(gs[1,:])
    levels = np.arange(0,1.51,0.01) 
    if day == '26':
        im.plt_ce_image(fig,ax1,time_MRR[:-1],height_MRR, SWC,levels, 
             v_min=0., v_max=1.5,xmax=time_MRR.max(), ymax=height_MRR.max(), 
             xticks=im.xticks1, yticks=im.yticks1,
             cb_ticks=levels[::30], var_label='SWC [g m$^{-3}$]')
    else:
        im.plt_ce_image(fig,ax1,time_MRR, height_MRR, SWC, levels, 
             v_min=0., v_max=1.5,xmax=time_MRR.max(), ymax=height_MRR.max(), 
             xticks=im.xticks1, yticks=im.yticks1,
             cb_ticks=levels[::30], var_label='SWC$_{RETRIEVAL}$ [g m$^{-3}$]')
################
## 3rd subplot --- SNOW WATER CONTENT ML
    ax2 = plt.subplot(gs[2,:])
    levels = np.arange(0,1.51,0.01) 
    im.plt_ce_image(fig,ax2,np.transpose(time[0]), height[0], (ice_amount_ml[0][:,:]),levels, 
            v_min=0., v_max=1.5, xmax=24., ymax=3000., 
            xticks =im.xticks2, yticks=im.yticks1,
            cb_ticks=levels[::30], var_label='SWC$_{MEPS}$(EM0) [g m$^{-2}$]')
################
## 4th subplot --- SNOW WATER CONTENT ML/THICKNESS
    ax3 = plt.subplot(gs[3,:])
    levels = np.arange(0,0.0751,0.001) 
    im.plt_ce_image(fig,ax3,np.transpose(time[0]), height[0], (ice_amount_ml_TH[0][:,:]),levels, 
            v_min=0., v_max=0.075, xmax=24., ymax=3000., 
            xticks =im.xticks2, yticks=im.yticks1,
            cb_ticks=levels[::15], var_label='SWC$_{MEPS}$(EM0) [g m$^{-3}$]')
################
## 5th subplot --- SWP
    ax4 = plt.subplot(gs[4,:])
# ensemble members
    for ens_memb in range(2,10):
        ax4.plot(time[ens_memb][0,:], SWP_model[ens_memb], color = im.memb_col, linestyle='-', label='_nolegend_')
    ax4.plot(time[1][0,:], SWP_model[1], color = im.memb_col, linestyle='-', label='ensemble member [g m$^{-2}$]')
    ax4.plot(time[0][0,:], SWP_model[0], 'k', linewidth=4, label='best guess [g m$^{-2}$]')
    
    if day == '26':
        line_ret = ax4.plot(time_MRR[:-1]/(60*60), SWP_ret, label = 'retrieval [g m$^{-2}$]')
    else:
        line_ret = ax4.plot(time_MRR/(60*60), SWP_ret, label = 'retrieval [g m$^{-2}$]')
    ax4.axis([0, 24, 0, 6000.])
### fine tuning
    lgd = plt.legend(loc='upper left',fontsize=im.label_fs)
    plt.grid()
# labels
    im.labels_x(ax4,im.xticks2)
    im.labels_y(ax4,im.yticks2,im.yl2, 'SWP')
    
    
    plt.tight_layout(pad=1.4, w_pad=1.5, h_pad=1.5)

    if savefig ==1:
        sF.save_figure_portrait(figdir, fig_name, form)
        print('saved: %s/%s' %(figdir, fig_name))
    else:
        plt.show()
    plt.close()
    
    for ens_memb in range(0,10):
        fn_snow[ens_memb].close()
        fn_graupel[ens_memb].close()
        
#    print(np.nanmax(ice_amount_ml_TH[0]))<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup