
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
import xarray as xr

import numpy as np
import metpy
import metpy.calc as mpcalc
from metpy.plots import ctables
from metpy.cbook import get_test_data
from metpy.units import units
import os
import scipy.integrate as integrate
import datetime as dt
import glob
import json
from datetime import datetime
from datetime import timedelta
from metpy.plots import ctables
from matplotlib.colors import Normalize
from matplotlib.colors import ListedColormap, LinearSegmentedColormap, BoundaryNorm

#import wrf
import scipy
#import xcape
import xarray


file_dir = '/data/icond2/'
os.chdir(file_dir)

data_det_ceil = xr.open_dataset('icond2_ceiling.nc').sel(longitude=slice(3,9,1),latitude=slice(49,51,1))

file_dir = '/data/icond2eps/hsurf'
os.chdir(file_dir)

data_hsurf = xr.open_dataset('icond2_hsurf.nc').sel(longitude=slice(3,9,1),latitude=slice(49,51,1))

file_dir = '/data/icond2eps/'
os.chdir(file_dir)

data = xarray.open_dataset('icond2eps_ceiling_latlon.grib2', engine='cfgrib')
data.to_netcdf('icond2eps_ceiling.nc')
data_ceil = xr.open_dataset('icond2eps_ceiling.nc').sel(longitude=slice(3,9,1),latitude=slice(49,51.02,1))

print(data_det_ceil)
print(data_ceil)
print(data_hsurf)

Ignoring index file 'icond2eps_ceiling_latlon.grib2.923a8.idx' older than GRIB file

<xarray.Dataset>
Dimensions:    (longitude: 301, latitude: 101, time: 48)
Coordinates:
  * longitude  (longitude) float32 3.0 3.02 3.04 3.06 ... 8.94 8.96 8.98 9.0
  * latitude   (latitude) float32 49.0 49.02 49.04 49.06 ... 50.96 50.98 51.0
  * time       (time) datetime64[ns] 2023-07-13T04:00:00 ... 2023-07-15T03:00:00
Data variables:
    ceil       (time, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2023-07-13 04:43:19 GMT by grib_to_netcdf-2.6.0: grib_to_ne...
<xarray.Dataset>
Dimensions:     (number: 20, step: 49, latitude: 101, longitude: 301)
Coordinates:
  * number      (number) int64 1 2 3 4 5 6 7 8 9 ... 12 13 14 15 16 17 18 19 20
    time        datetime64[ns] ...
  * step        (step) timedelta64[ns] 00:00:00 01:00:00 ... 2 days 00:00:00
    level       float64 ...
  * latitude    (latitude) float64 49.0 49.02 49.04 49.06 ... 50.96 50.98 51.0
  * longitude   (longitude) float64 3.0 3.02 3.04 3.06 ... 8.94 8.96 8.98 9.0
    valid_time  (step) datetime64[ns] ...
Data variables:
    ceil        (number, step, latitude, longitude) float32 ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             edzw
    GRIB_centreDescription:  Offenbach 
    GRIB_subCentre:          255
    Conventions:             CF-1.7
    institution:             Offenbach 
    history:                 2023-07-13T06:03 GRIB to CDM+CF via cfgrib-0.9.9...
<xarray.Dataset>
Dimensions:    (longitude: 301, latitude: 101, time: 1)
Coordinates:
  * longitude  (longitude) float32 3.0 3.02 3.04 3.06 ... 8.94 8.96 8.98 9.0
  * latitude   (latitude) float32 49.0 49.02 49.04 49.06 ... 50.96 50.98 51.0
  * time       (time) datetime64[ns] 2021-12-22T09:00:00
Data variables:
    h          (time, latitude, longitude) float32 ...
Attributes:
    Conventions:  CF-1.6
    history:      2021-12-22 19:08:59 GMT by grib_to_netcdf-2.6.0: grib_to_ne...


# To parse the full dataset, we can call parse_cf without an argument, and assign the returned Dataset.

data_ceil = data_ceil.metpy.parse_cf()
data_det_ceil = data_det_ceil.metpy.parse_cf()

x, y = data_ceil['ceil'].metpy.coordinates('x', 'y')

time = data_ceil['ceil'].step
member = data_ceil['ceil'].number
#time2 = data_det_rr['tp'].metpy.time

timeinit = data_ceil.time
timeinit = datetime.utcfromtimestamp(timeinit.item()/1e9)
print(timeinit)

hsurf = data_hsurf['h']*3.28084 #convert to feet
#uh = data_uh['UH_MAX']
ceil = data_ceil['ceil']
#rr = data_rr['tp']
#det_rr = data_det_rr['tp']
det_ceil = data_det_ceil['ceil']
#det_uh = data_det_uhmax['UH_MAX']

#rr.data = np.nan_to_num(rr.data, copy=True, nan=0)
print(np.shape(ceil.data))

#vmax_median = np.percentile(vmax, 50)
ceil_median = np.empty((49,101,301))
ceil_95 = np.empty((49,101,301))

for i in range(0,49):
    for j in range (0,101):
        for k in range(0,301):
            ceil_median[i,j,k] = np.percentile(ceil.data[:,i,j,k],50)

print(np.shape(ceil_median))

for i in range(0,49):
    for j in range (0,101):
        for k in range(0,301):
            ceil_95[i,j,k] = np.percentile(ceil.data[:,i,j,k],5)

/home/lmathias/anaconda3/envs/metpy/lib/python3.9/site-packages/metpy/xarray.py:349: UserWarning: More than one time coordinate present for variable "ceil".
  warnings.warn('More than one ' + axis + ' coordinate present for variable'
Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable
Found valid latitude/longitude coordinates, assuming latitude_longitude for projection grid_mapping variable

2023-07-13 03:00:00
(20, 49, 101, 301)
(49, 101, 301)


ceil_prob_15 = np.empty((49,101,301))
ceil_prob_10 = np.empty((49,101,301))
ceil_prob_5 = np.empty((49,101,301))
ceil_prob_2 = np.empty((49,101,301))
ceil_agl = np.empty((20,49,101,301))

z = 1500
v = 1000
w = 500
n = 200

ceil = data_ceil['ceil']*3.28084

for i in range(0,20):
    for j in range(0,49):
        for k in range (0,101):
            for l in range(0,301):
                ceil_agl[i,j,k,l] = ceil.data[i,j,k,l]-hsurf.data[0,k,l]


# for i in range(0,49):
#     for j in range (0,101):
#         for k in range(0,301):
#             ceil_prob_15[i,j,k] = (sum(h < z for h in ceil_agl[:,i,j,k])/20)*100
            
# for i in range(0,49):
#     for j in range (0,101):
#         for k in range(0,301):
#             ceil_prob_10[i,j,k] = (sum(h < v for h in ceil_agl[:,i,j,k])/20)*100
            
# for i in range(0,49):
#     for j in range (0,101):
#         for k in range(0,301):
#             ceil_prob_5[i,j,k] = (sum(h < w for h in ceil_agl[:,i,j,k])/20)*100
            
# for i in range(0,49):
#     for j in range (0,101):
#         for k in range(0,301):
#             ceil_prob_2[i,j,k] = (sum(h < n for h in ceil_agl[:,i,j,k])/20)*100

            
for i in range(0,49):
    for j in range (0,101):
        for k in range(0,301):
            count = np.count_nonzero((ceil_agl[:,i,j,k]) < 1500)
            ceil_prob_15[i,j,k] = (count/20)*100
            count = 0
            
for i in range(0,49):
    for j in range (0,101):
        for k in range(0,301):
            count = np.count_nonzero((ceil_agl[:,i,j,k]) < 1000)
            ceil_prob_10[i,j,k] = (count/20)*100
            count = 0
            
for i in range(0,49):
    for j in range (0,101):
        for k in range(0,301):
            count = np.count_nonzero((ceil_agl[:,i,j,k]) < 500)
            ceil_prob_5[i,j,k] = (count/20)*100
            count = 0
            
for i in range(0,49):
    for j in range (0,101):
        for k in range(0,301):
            count = np.count_nonzero((ceil_agl[:,i,j,k]) < 200)
            ceil_prob_2[i,j,k] = (count/20)*100
            count = 0

#print(ceil_prob_2)


def plot_background(ax):
    ax.set_extent([5, 7, 49.1, 50.5])
    ax.add_feature(cfeature.COASTLINE.with_scale('10m'), LineWidth=2)
    ax.add_feature(cfeature.BORDERS.with_scale('10m'),LineWidth=2)
    #gl = ax.gridlines(draw_labels=True,linewidth=0.5, color='gray', alpha=0.5, linestyle='--')
    #gl.xlabels_top = False
    #gl.ylabels_right = False
    #gl.xlabel_style = {'size': 12, 'color': 'black'}
    #gl.ylabel_style = {'size': 12, 'color': 'black'}
    #gl.xformatter = LONGITUDE_FORMATTER
    #gl.yformatter = LATITUDE_FORMATTER
    return ax

#import matplotlib
#cmap = matplotlib.cm.get_cmap('cubehelix_r')
#for i in range(20):
    #rgba = cmap(i)
    # rgb2hex accepts rgb or rgba
    #print(rgba)
    
cmap = ctables.colortables.get_colortable('NWSStormClearReflectivity')
newcmap = ListedColormap(cmap(np.linspace(0.25, 0.92, 28)))

cmap2 = ctables.colortables.get_colortable('NWSReflectivity')
newcmap2 = ListedColormap(cmap2(np.linspace(0.2, 0.96, 28)))

cmap4 = ctables.colortables.get_colortable('precipitation')
newcmap4 = ListedColormap(cmap4(np.linspace(0, 0.75, 15)))


bounds = [0,0.1,0.5,1,3,5,10,15,20,25,30,35,40,45,50,55]
norm = BoundaryNorm(bounds, newcmap4.N)
bounds6 = [0,0.1,5,10,15,20,25,30,40,50,60,70,80,90,100,110]
norm6 = BoundaryNorm(bounds6, newcmap4.N)
bounds24 = [0,0.1,5,10,20,30,40,50,60,70,80,90,100,110,120,130]
norm24 = BoundaryNorm(bounds6, newcmap4.N)

cmap5 = plt.cm.get_cmap('Dark2_r')
#newcmap5 = ListedColormap(cmap5(np.linspace(0, 1, 8)))
bounds_ceil = [0,100,200,300,400,500,750,1000,1500]
norm_ceil = BoundaryNorm(bounds_ceil, cmap5.N)

#cmap2 = ListedColormap(colors2)
#newcmap2 = ListedColormap(cmap2(np.linspace(0, 0.9, 29)))


# Create the figure and plot background on different axes
crs = ccrs.Mercator()

for i in range(1,49):
    
    fig, axarr = plt.subplots(nrows=1, ncols=3, figsize=(25, 10), constrained_layout=False,
                  subplot_kw={'projection': crs})
    # Set height padding for plots
    fig.set_constrained_layout_pads(w_pad=0., h_pad=10, hspace=0., wspace=0.)
    
    axlist = axarr.flatten()
    
    for ax in axlist:
        plot_background(ax)
        
    timestep=timeinit+timedelta(hours=i)
    time2 = data_det_ceil['ceil'].metpy.time
    time3 = data_hsurf['h'].metpy.time
    clevs_ceil= np.arange(0,1500,100)
    # cmap = plt.get_cmap('gist_ncar')
    # newcmap = ListedColormap(cmap(np.linspace(0.15, 0.9, 30)))
        
    # Upper left plot
    cf1 = axlist[0].contourf(data_ceil.longitude, data_ceil.latitude, (ceil_median[i,:,:]*3.28084)-hsurf.metpy.loc[{'time': time3[0]}], 
                      [0,100,200,300,400,500,750,1000,1500], cmap=cmap5, norm=norm_ceil,transform=ccrs.PlateCarree())
    #ccf1= axlist[0].contour(data_ceil.longitude, data_ceil.latitude,  ceil_median[i,:,:],
                     #[100,200,500,1000,1500], colors='dimgrey', linestyles="dotted",transform=ccrs.PlateCarree())
    #axlist[0].clabel(ccf1, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
    axlist[0].set_title('50th Percentile (Median)', fontsize=16)
    #cb1= fig.colorbar(cf1, ax=axlist[0], orientation='vertical', 
                     # ticks=(10,20,30,40,50,60,70,80,90,100,110,120,130),
                            #  shrink=0.73, fraction=0.1, pad=0)
    #cb1.set_label('km/h', size='x-large')
    
    cf2 = axlist[1].contourf(data_ceil.longitude, data_ceil.latitude, (ceil_95[i,:,:]*3.28084)-hsurf.metpy.loc[{'time': time3[0]}], 
                      [0,100,200,300,400,500,750,1000,1500], cmap=cmap5, norm=norm_ceil,transform=ccrs.PlateCarree())
    #ccf2= axlist[1].contour(data_ceil.longitude, data_ceil.latitude,  ceil_95[i,:,:],
                     #[100,200,500,1000,1500], colors='dimgrey', linestyles="dotted",transform=ccrs.PlateCarree())
    #axlist[1].clabel(ccf2, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
    axlist[1].set_title('5th Percentile', fontsize=16)
    #cb2= fig.colorbar(cf2, ax=axlist[1], orientation='vertical', 
                      #ticks=(10,20,30,40,50,60,70,80,90,100,110,120,130),
                             # shrink=0.73, fraction=0.1, pad=0)
   # cb2.set_label('km/h', size='x-large')

    cf3 = axlist[2].contourf(data_det_ceil.longitude, data_det_ceil.latitude, (det_ceil.metpy.loc[{'time': time2[i-1]}]*3.28084)-hsurf.metpy.loc[{'time': time3[0]}], 
                          [0,100,200,300,400,500,750,1000,1500],cmap=cmap5, norm=norm_ceil,transform=ccrs.PlateCarree())
    #ccf3= axlist[2].contour(data_det_vmax.longitude, data_det_vmax.latitude, det_gust.metpy.loc[{'time': time2[i-1]}]*3.6,
                     #[30,50,70,90,110,130], colors='dimgrey', linestyles="dotted",transform=ccrs.PlateCarree())
    #axlist[2].clabel(ccf3, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
    axlist[2].set_title('Deterministic', fontsize=16)
   # cb3= fig.colorbar(cf3, ax=axlist[2], orientation='vertical', 
                      #ticks=(10,20,30,40,50,60,70,80,90,100,110,120,130),
                              #shrink=0.73, fraction=0.1, pad=0)
    #cb3.set_label('km/h', size='x-large')
    cb = fig.colorbar(cf1, ax=axarr.ravel().tolist(), orientation='vertical',
                      ticks=(0,100,200,300,400,500,750,1000,1500), fraction=0.01, aspect=30, pad=0.02)
    cb.set_label('ft AGL', size='x-large')


    # Set figure title
    plt.gcf().text(0.130, 0.90, 'Model: ICON-D2-EPS 0.02° | ' + timeinit.strftime('Init: %d.%m.%Y %H:%M UTC | ')+timestep.strftime('Valid: %d.%m.%Y %H:%M UTC'), fontsize=20)
    plt.gcf().text(0.130, 0.86, 'Parameter: Ceiling', fontsize=20)
    # Display the plot
    time2 = str(i*1)
    base_filename='icond2eps_ceil_perc_'
    suffix='.jpeg'
    my_file = base_filename+time2+suffix
    print(my_file)
    plt.savefig(my_file, format="jpeg", bbox_inches='tight', dpi=85)
    plt.close(fig)

icond2eps_ceil_perc_1.jpeg

/home/lmathias/anaconda3/envs/metpy/lib/python3.9/site-packages/cartopy/mpl/feature_artist.py:211: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later
  c = matplotlib.collections.PathCollection(paths,

icond2eps_ceil_perc_2.jpeg

/home/lmathias/anaconda3/envs/metpy/lib/python3.9/site-packages/cartopy/mpl/feature_artist.py:211: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later
  c = matplotlib.collections.PathCollection(paths,

icond2eps_ceil_perc_3.jpeg

/home/lmathias/anaconda3/envs/metpy/lib/python3.9/site-packages/cartopy/mpl/feature_artist.py:211: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later
  c = matplotlib.collections.PathCollection(paths,

icond2eps_ceil_perc_4.jpeg

/home/lmathias/anaconda3/envs/metpy/lib/python3.9/site-packages/cartopy/mpl/feature_artist.py:211: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later
  c = matplotlib.collections.PathCollection(paths,

icond2eps_ceil_perc_5.jpeg

/home/lmathias/anaconda3/envs/metpy/lib/python3.9/site-packages/cartopy/mpl/feature_artist.py:211: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later
  c = matplotlib.collections.PathCollection(paths,


# Create the figure and plot background on different axes
crs = ccrs.Mercator()
ceil_prob_15_smooth = scipy.ndimage.zoom(ceil_prob_15, 1.2)
ceil_prob_10_smooth = scipy.ndimage.zoom(ceil_prob_10, 1.2)
ceil_prob_5_smooth = scipy.ndimage.zoom(ceil_prob_5, 1.2)
ceil_prob_2_smooth = scipy.ndimage.zoom(ceil_prob_2, 1.2)
x_smooth = scipy.ndimage.zoom(x, 1.2)
y_smooth = scipy.ndimage.zoom(y, 1.2)
ceil_prob_15_smooth = np.where(ceil_prob_15_smooth < 100, ceil_prob_15_smooth, 100)
ceil_prob_10_smooth = np.where(ceil_prob_10_smooth < 100, ceil_prob_10_smooth, 100)
ceil_prob_5_smooth = np.where(ceil_prob_5_smooth < 100, ceil_prob_5_smooth, 100)
ceil_prob_2_smooth = np.where(ceil_prob_2_smooth < 100, ceil_prob_2_smooth, 100)
ceil_prob_15_smooth = np.where(ceil_prob_15_smooth > 0, ceil_prob_15_smooth, 0)
ceil_prob_10_smooth = np.where(ceil_prob_10_smooth > 0, ceil_prob_10_smooth, 0)
ceil_prob_5_smooth = np.where(ceil_prob_5_smooth > 0, ceil_prob_5_smooth, 0)
ceil_prob_2_smooth = np.where(ceil_prob_2_smooth > 0, ceil_prob_2_smooth, 0)

#print(np.max(ceil_prob_15_smooth))

for i in range(1,49):
    
    fig, axarr = plt.subplots(nrows=1, ncols=4, figsize=(25, 10), constrained_layout=False,
                  subplot_kw={'projection': crs})
    # Set height padding for plots
    fig.set_constrained_layout_pads(w_pad=0., h_pad=10, hspace=0., wspace=0.)
    
    axlist = axarr.flatten()
    
    for ax in axlist:
        plot_background(ax)
        
    timestep=timeinit+timedelta(hours=i)
    time2 = data_det_ceil['ceil'].metpy.time
    time3 = data_hsurf['h'].metpy.time
    clevs_prob= np.arange(0,110,10)
    # cmap = plt.get_cmap('gist_ncar')
    # newcmap = ListedColormap(cmap(np.linspace(0.15, 0.9, 30)))
        
    # Upper left plot
    cf1 = axlist[0].contourf(x, y, ceil_prob_15[i,:,:],
                      clevs_prob, cmap='BuPu', transform=ccrs.PlateCarree())
    #ccf1= axlist[0].contour(data_ceil.longitude, data_ceil.latitude,  ceil_median[i,:,:],
                     #[100,200,500,1000,1500], colors='dimgrey', linestyles="dotted",transform=ccrs.PlateCarree())
    #axlist[0].clabel(ccf1, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
    axlist[0].set_title('Probability h < 1500 ft AGL', fontsize=16)
    #cb1= fig.colorbar(cf1, ax=axlist[0], orientation='vertical', 
                     # ticks=(10,20,30,40,50,60,70,80,90,100,110,120,130),
                            #  shrink=0.73, fraction=0.1, pad=0)
    #cb1.set_label('km/h', size='x-large')
    
    cf2 = axlist[1].contourf(x, y, ceil_prob_10[i,:,:],
                      clevs_prob, cmap='BuPu', transform=ccrs.PlateCarree())
    #ccf2= axlist[1].contour(data_ceil.longitude, data_ceil.latitude,  ceil_95[i,:,:],
                     #[100,200,500,1000,1500], colors='dimgrey', linestyles="dotted",transform=ccrs.PlateCarree())
    #axlist[1].clabel(ccf2, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
    axlist[1].set_title('Probability h < 1000 ft AGL', fontsize=16)
    #cb2= fig.colorbar(cf2, ax=axlist[1], orientation='vertical', 
                      #ticks=(10,20,30,40,50,60,70,80,90,100,110,120,130),
                             # shrink=0.73, fraction=0.1, pad=0)
   # cb2.set_label('km/h', size='x-large')

    cf3 = axlist[2].contourf(x, y, ceil_prob_5[i,:,:], 
                      clevs_prob, cmap='BuPu', transform=ccrs.PlateCarree())
    #ccf3= axlist[2].contour(data_det_vmax.longitude, data_det_vmax.latitude, det_gust.metpy.loc[{'time': time2[i-1]}]*3.6,
                     #[30,50,70,90,110,130], colors='dimgrey', linestyles="dotted",transform=ccrs.PlateCarree())
    #axlist[2].clabel(ccf3, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
    axlist[2].set_title('Probability h < 500 ft AGL', fontsize=16)
   # cb3= fig.colorbar(cf3, ax=axlist[2], orientation='vertical', 
                      #ticks=(10,20,30,40,50,60,70,80,90,100,110,120,130),
                              #shrink=0.73, fraction=0.1, pad=0)
    #cb3.set_label('km/h', size='x-large')

    cf4 = axlist[3].contourf(x, y, ceil_prob_2[i,:,:], 
                      clevs_prob, cmap='BuPu',transform=ccrs.PlateCarree())
    #ccf3= axlist[2].contour(data_det_vmax.longitude, data_det_vmax.latitude, det_gust.metpy.loc[{'time': time2[i-1]}]*3.6,
                     #[30,50,70,90,110,130], colors='dimgrey', linestyles="dotted",transform=ccrs.PlateCarree())
    #axlist[2].clabel(ccf3, fontsize=10, inline=1, inline_spacing=1, fmt='%i', rightside_up=True)
    axlist[3].set_title('Probability h < 200 ft AGL', fontsize=16)
    cb = fig.colorbar(cf1, ax=axarr.ravel().tolist(), orientation='vertical',
                      ticks=(0,10,20,30,40,50,60,70,80,90,100), fraction=0.0075, aspect=30, pad=0.02)
    cb.set_label('%', size='x-large')


    # Set figure title
    plt.gcf().text(0.130, 0.86, 'Model: ICON-D2-EPS 0.02° | ' + timeinit.strftime('Init: %d.%m.%Y %H:%M UTC | ')+timestep.strftime('Valid: %d.%m.%Y %H:%M UTC'), fontsize=20)
    plt.gcf().text(0.130, 0.82, 'Parameter: Ceiling', fontsize=20)
    # Display the plot
    time2 = str(i*1)
    base_filename='icond2eps_ceil_prob_'
    suffix='.jpeg'
    my_file = base_filename+time2+suffix
    print(my_file)
    plt.savefig(my_file, format="jpeg", bbox_inches='tight', dpi=85)
    plt.close(fig)

icond2eps_ceil_prob_1.jpeg
icond2eps_ceil_prob_2.jpeg
icond2eps_ceil_prob_3.jpeg
icond2eps_ceil_prob_4.jpeg
icond2eps_ceil_prob_5.jpeg
icond2eps_ceil_prob_6.jpeg
icond2eps_ceil_prob_7.jpeg
icond2eps_ceil_prob_8.jpeg
icond2eps_ceil_prob_9.jpeg
icond2eps_ceil_prob_10.jpeg
icond2eps_ceil_prob_11.jpeg
icond2eps_ceil_prob_12.jpeg
icond2eps_ceil_prob_13.jpeg
icond2eps_ceil_prob_14.jpeg
icond2eps_ceil_prob_15.jpeg
icond2eps_ceil_prob_16.jpeg
icond2eps_ceil_prob_17.jpeg
icond2eps_ceil_prob_18.jpeg
icond2eps_ceil_prob_19.jpeg
icond2eps_ceil_prob_20.jpeg
icond2eps_ceil_prob_21.jpeg
icond2eps_ceil_prob_22.jpeg
icond2eps_ceil_prob_23.jpeg
icond2eps_ceil_prob_24.jpeg
icond2eps_ceil_prob_25.jpeg
icond2eps_ceil_prob_26.jpeg
icond2eps_ceil_prob_27.jpeg
icond2eps_ceil_prob_28.jpeg
icond2eps_ceil_prob_29.jpeg
icond2eps_ceil_prob_30.jpeg
icond2eps_ceil_prob_31.jpeg
icond2eps_ceil_prob_32.jpeg
icond2eps_ceil_prob_33.jpeg
icond2eps_ceil_prob_34.jpeg
icond2eps_ceil_prob_35.jpeg
icond2eps_ceil_prob_36.jpeg
icond2eps_ceil_prob_37.jpeg
icond2eps_ceil_prob_38.jpeg
icond2eps_ceil_prob_39.jpeg
icond2eps_ceil_prob_40.jpeg
icond2eps_ceil_prob_41.jpeg
icond2eps_ceil_prob_42.jpeg
icond2eps_ceil_prob_43.jpeg
icond2eps_ceil_prob_44.jpeg
icond2eps_ceil_prob_45.jpeg
icond2eps_ceil_prob_46.jpeg
icond2eps_ceil_prob_47.jpeg
icond2eps_ceil_prob_48.jpeg