Notebook Source Codes

Share your source codes

Source codes: 125

Un exemple 1D de code python sur la méthode des k plus proches voisins

					
import matplotlib.pyplot as plt
import random 
import math

#----------------------------------------------------------------------------------------#
# Create a dataset

x_list = [0.001*random.randrange(0, 10000) for i in range(20)]
y_list = [1 if x > 5 else -1 for x in x_list]
y_target_list = [-1,1]

#----------------------------------------------------------------------------------------#
# DataSet Viz.

plt.scatter(x_list,y_list)

plt.xlabel('x')
plt.ylabel('y target')

plt.title('k-Nearest Neighbors 1D')

plt.savefig("k_nearest_neighbors_1d.png", bbox_inches='tight')
plt.show()

#----------------------------------------------------------------------------------------#
# For a new point x_new; do a prediction y

def prediction(x_new,k):
	# calculate the distance for each point in the dataset
	point_list = []
	for idx,x in enumerate(x_list):
		d = math.fabs(x - x_new)
		y = y_list[idx]
		point_list.append((d,y))
	# sort point_list by distance to x_new
	point_list = sorted(point_list, key=lambda x: x[0])
	# keep only k nearest point(s)
	if k < len(point_list): 
		k_nearest_points_list = point_list[:k]
	else:	
		k_nearest_points_list = point_list
	# prediction step 1
	y_new_count_list = [0 for y in y_target_list]
	for k_nearest_point in k_nearest_points_list:
		y_idx = y_target_list.index(k_nearest_point[1])
		y_new_count_list[y_idx] = y_new_count_list[y_idx] + 1
	max_count = max(y_new_count_list)
	idx_max_list = [idx for idx,y_new_count in enumerate(y_new_count_list) if y_new_count == max_count]
	# prediction step 2
	if len(idx_max_list) == 1:
		y_new = y_target_list[idx_max_list[0]]
	else:
		y_new = [y_target_list[idx] for idx in idx_max_list]	
	# return
	return y_new

#----------------------------------------------------------------------------------------#
# Test 1

k = 1
x_new = 3
y_new = prediction(x_new,k)

print('Test 1: ',x_new,y_new)

#----------------------------------------------------------------------------------------#
# Test 2

k = 1
x_new = 7
y_new = prediction(x_new,k)

print('Test 2: ',x_new,y_new)

#----------------------------------------------------------------------------------------#
# Test 3

k = 2
x_new = 5
y_new = prediction(x_new,k)

print('Test 3: ',x_new,y_new)

#----------------------------------------------------------------------------------------#
# Test 4

k = 3
x_new = 7
y_new = prediction(x_new,k)

print('Test 4: ',x_new,y_new)

A k-nearest neighbors python algorithm for 1D cases

					
import matplotlib.pyplot as plt
import random 
import math

#----------------------------------------------------------------------------------------#
# Create a dataset

x_list = [0.001*random.randrange(0, 10000) for i in range(20)]
y_list = [1 if x > 5 else -1 for x in x_list]
y_target_list = [-1,1]

#----------------------------------------------------------------------------------------#
# DataSet Viz.

plt.scatter(x_list,y_list)

plt.xlabel('x')
plt.ylabel('y target')

plt.title('k-Nearest Neighbors 1D')

plt.savefig("k_nearest_neighbors_1d.png", bbox_inches='tight')
plt.show()

#----------------------------------------------------------------------------------------#
# For a new point x_new; do a prediction y

def prediction(x_new,k):
	# calculate the distance for each point in the dataset
	point_list = []
	for idx,x in enumerate(x_list):
		d = math.fabs(x - x_new)
		y = y_list[idx]
		point_list.append((d,y))
	# sort point_list by distance to x_new
	point_list = sorted(point_list, key=lambda x: x[0])
	# keep only k nearest point(s)
	if k < len(point_list): 
		k_nearest_points_list = point_list[:k]
	else:	
		k_nearest_points_list = point_list
	# prediction step 1
	y_new_count_list = [0 for y in y_target_list]
	for k_nearest_point in k_nearest_points_list:
		y_idx = y_target_list.index(k_nearest_point[1])
		y_new_count_list[y_idx] = y_new_count_list[y_idx] + 1
	max_count = max(y_new_count_list)
	idx_max_list = [idx for idx,y_new_count in enumerate(y_new_count_list) if y_new_count == max_count]
	# prediction step 2
	if len(idx_max_list) == 1:
		y_new = y_target_list[idx_max_list[0]]
	else:
		y_new = [y_target_list[idx] for idx in idx_max_list]	
	# return
	return y_new

#----------------------------------------------------------------------------------------#
# Test 1

k = 1
x_new = 3
y_new = prediction(x_new,k)

print('Test 1: ',x_new,y_new)

#----------------------------------------------------------------------------------------#
# Test 2

k = 1
x_new = 7
y_new = prediction(x_new,k)

print('Test 2: ',x_new,y_new)

#----------------------------------------------------------------------------------------#
# Test 3

k = 2
x_new = 5
y_new = prediction(x_new,k)

print('Test 3: ',x_new,y_new)

#----------------------------------------------------------------------------------------#
# Test 4

k = 3
x_new = 7
y_new = prediction(x_new,k)

print('Test 4: ',x_new,y_new)

Plot a year of MODIS L2 Gridded Multilayer Clouds (all pixels) with python 3

					
#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm
from pyhdf.SD import SD, SDC 

import pprint
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import calendar 

year = 2008

file_name = 'gridded_data_{:4d}_{:02d}.txt'.format(year,1)

date = str(year)

input_data = np.loadtxt(file_name, skiprows=1)

input_data_shape = input_data.shape

dict_settings = {}

#----------------------------------------------------------------------------------------#
# Full Year

input_data = np.zeros((input_data_shape))

for month in range(1,13):
	file_name = 'gridded_data_{:4d}_{:02d}.txt'.format(year,month)
	input_data_month = np.loadtxt(file_name, skiprows=1)
	for i in range(input_data_shape[0]):
		for j in range(input_data_shape[1]):
			input_data[i,j] = input_data[i,j] + input_data_month[i,j]
	
#----------------------------------------------------------------------------------------#
# Plot gidded data on a global map function

def plot_gridded_data(data,dict_settings):
	
	fig = plt.figure()
	
	m = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90, 
				llcrnrlon=-180,urcrnrlon=180,resolution='c')
	
	m.drawcoastlines()
	m.drawparallels(np.arange(-90.,90.,30.))
	m.drawmeridians(np.arange(-180.,180.,30.))
	
	plt.title(dict_settings['fig_title'], size=8)
	
	cmap = [(0.75,0.75,0.75)] + [(cm.jet (i)) for i in range(1,256)] 
	cmap = mpl.colors.ListedColormap(cmap)
	
	bounds = dict_settings['bounds']
	
	norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
	
	img = m.imshow(data, cmap=cmap, norm=norm, interpolation='none')
	
	ticks = dict_settings['ticks']
		
	cb = m.colorbar(img,"bottom", size="5%", pad='4%', boundaries=bounds, 
					ticks=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])
	
	cb.ax.tick_params(labelsize=6) 
	
	plt.savefig(dict_settings['output_file_name'], bbox_inches='tight', dpi=200)
	plt.close()

#----------------------------------------------------------------------------------------#
# plot cloud fraction

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,1] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Cloudy\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_global_map_modis_l2_gridded_data_cloud_mask.png'.format(year)
dict_settings['bounds'] = [-1.0,0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]
dict_settings['ticks'] = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud 

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,2] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_global_map_modis_l2_gridded_data_c6_multilayer_clouds.png'.format(year)
dict_settings['bounds'] = [-1.0,0.0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.40,0.45,0.5]
dict_settings['ticks'] = [0, 0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.40,0.45,0.5]

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 (without PH) multilayer cloud 

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,3] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ without\ PH\ test\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_global_map_modis_l2_gridded_data_c6_multilayer_clouds_without_PH.png'.format(year)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C5 multilayer cloud 

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,4] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MMODIS C5 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_global_map_modis_l2_gridded_data_c5_multilayer_clouds.png'.format(year)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud QA 1

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,5] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ -\ QA1\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_global_map_modis_l2_gridded_data_c6_multilayer_clouds_qa1.png'.format(year)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud QA 2

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,6] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ -\ QA2\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_global_map_modis_l2_gridded_data_c6_multilayer_clouds_qa2.png'.format(year)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud QA 3

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,7] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ -\ QA3\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_global_map_modis_l2_gridded_data_c6_multilayer_clouds_qa3.png'.format(year)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud QA 4

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,8] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ -\ QA4\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_global_map_modis_l2_gridded_data_c6_multilayer_clouds_qa4.png'.format(year)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud QA 5

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,9] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ -\ QA5\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_global_map_modis_l2_gridded_data_c6_multilayer_clouds_qa5.png'.format(year)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#

Plot a month of MODIS L2 Gridded Multilayer Clouds (all pixels) with python 3

					
#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm
from pyhdf.SD import SD, SDC 

import pprint
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import calendar 

year = 2008
month = 12

file_name = 'gridded_data_{:4d}_{:02d}.txt'.format(year,month)

date = calendar.month_name[month] + ' ' + str(year)

input_data = np.loadtxt(file_name, skiprows=1)

print(input_data.shape)

dict_settings = {}

#----------------------------------------------------------------------------------------#
# Plot gidded data on a global map function

def plot_gridded_data(data,dict_settings):
	
	fig = plt.figure()
	
	m = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90, 
				llcrnrlon=-180,urcrnrlon=180,resolution='c')
	
	m.drawcoastlines()
	m.drawparallels(np.arange(-90.,90.,30.))
	m.drawmeridians(np.arange(-180.,180.,30.))
	
	plt.title(dict_settings['fig_title'], size=8)
	
	cmap = [(0.75,0.75,0.75)] + [(cm.jet (i)) for i in range(1,256)] 
	cmap = mpl.colors.ListedColormap(cmap)
	
	bounds = [-1.0,0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]
	
	norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
	
	img = m.imshow(data, cmap=cmap, norm=norm, interpolation='none')
	
	bounds = [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]
	
	cb = m.colorbar(img,"bottom", size="5%", pad='4%', boundaries=bounds, 
					ticks=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])
	
	cb.ax.tick_params(labelsize=6) 
	
	plt.savefig(dict_settings['output_file_name'], bbox_inches='tight', dpi=200)
	plt.close()

#----------------------------------------------------------------------------------------#
# plot cloud fraction

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,1] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Cloudy\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_{:02d}_global_map_monthly_modis_l2_gridded_data_cloud_mask.png'.format(year,month)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud 

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,2] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_{:02d}_global_map_monthly_modis_l2_gridded_data_c6_multilayer_clouds.png'.format(year,month)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 (without PH) multilayer cloud 

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,3] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ without\ PH\ test\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_{:02d}_global_map_monthly_modis_l2_gridded_data_c6_multilayer_clouds_without_PH.png'.format(year,month)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C5 multilayer cloud 

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,4] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MMODIS C5 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_{:02d}_global_map_monthly_modis_l2_gridded_data_c5_multilayer_clouds.png'.format(year,month)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud QA 1

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,5] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ -\ QA1\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_{:02d}_global_map_monthly_modis_l2_gridded_data_c6_multilayer_clouds_qa1.png'.format(year,month)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud QA 2

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,6] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ -\ QA2\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_{:02d}_global_map_monthly_modis_l2_gridded_data_c6_multilayer_clouds_qa2.png'.format(year,month)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud QA 3

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,7] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ -\ QA3\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_{:02d}_global_map_monthly_modis_l2_gridded_data_c6_multilayer_clouds_qa3.png'.format(year,month)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud QA 4

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,8] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ -\ QA4\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_{:02d}_global_map_monthly_modis_l2_gridded_data_c6_multilayer_clouds_qa4.png'.format(year,month)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#
# plot C6 multilayer cloud QA 5

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,9] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

dict_settings['fig_title'] = r'$P(\ Multilayer\ Cloud\ -\ QA5\ |\ 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - ' + date
dict_settings['output_file_name'] = '{:4d}_{:02d}_global_map_monthly_modis_l2_gridded_data_c6_multilayer_clouds_qa5.png'.format(year,month)

plot_gridded_data(data,dict_settings)

#----------------------------------------------------------------------------------------#

Plot MODIS L2 Gridded Data (all cloudy pixels) with python 3

					
#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm
from pyhdf.SD import SD, SDC 

import pprint
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm

file_name = 'gridded_data_2008_01.txt'

input_data = np.loadtxt(file_name, skiprows=1)

#----------------------------------------------------------------------------------------#
# plot cloud fraction

data = np.zeros((64800))

for i in range(64800):
	if input_data[i,0] > 0 and input_data[i,1] > 0:
		data[i] = input_data[i,1] /  input_data[i,0]
	else:
		data[i] = -0.5

data = data.reshape(180,360)

fig = plt.figure()

m = Basemap(projection='cyl',llcrnrlat=-90,urcrnrlat=90, 
			llcrnrlon=-180,urcrnrlon=180,resolution='c')

m.drawcoastlines()
m.drawparallels(np.arange(-90.,90.,30.))
m.drawmeridians(np.arange(-180.,180.,30.))

plt.title(r'$P(Cloudy | 1^o \times 1^o bin; Obs.)$'+' \n MODIS C6 Gridded Data L2 (all pixels) - January 2008', size=8)

cmap = [(0.75,0.75,0.75)] + [(cm.jet (i)) for i in range(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [-1.0,0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = m.imshow(data, cmap=cmap, norm=norm, interpolation='none')

bounds = [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]

cb = m.colorbar(img,"bottom", size="5%", pad='4%', boundaries=bounds, 
				ticks=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])

cb.ax.tick_params(labelsize=8) 

plt.savefig('global_map_monthly_modis_l2_gridded_data_cloud_mask.png', bbox_inches='tight', dpi=200)
plt.close()

#----------------------------------------------------------------------------------------#

How to plot a confusion matrix with matplotlib and seaborn

					
#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sn
import pandas as pd

import seaborn as sns
import math

from mpl_toolkits.axes_grid1 import make_axes_locatable

import matplotlib as mpl

mpl.style.use('seaborn')

conf_arr = np.array([[88,14,4],[12,85,11],[5,15,91]])

sum = conf_arr.sum()

conf_arr = conf_arr * 100.0 / ( 1.0 * sum )

df_cm = pd.DataFrame(conf_arr, 
  index = [ 'Dog', 'Cat', 'Rabbit'],
  columns = ['Dog', 'Cat', 'Rabbit'])

fig = plt.figure()

plt.clf()

ax = fig.add_subplot(111)
ax.set_aspect(1)

cmap = sns.cubehelix_palette(light=1, as_cmap=True)

res = sn.heatmap(df_cm, annot=True, vmin=0.0, vmax=100.0, fmt='.2f', cmap=cmap)

res.invert_yaxis()

plt.yticks([0.5,1.5,2.5], [ 'Dog', 'Cat', 'Rabbit'],va='center')

plt.title('Confusion Matrix')

plt.savefig('confusion_matrix.png', dpi=100, bbox_inches='tight' )

plt.close()

Tracer une matrice de confusion avec matplotlib et seaborn

					
#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sn
import pandas as pd

import seaborn as sns
import math

from mpl_toolkits.axes_grid1 import make_axes_locatable

import matplotlib as mpl

mpl.style.use('seaborn')

conf_arr = np.array([[88,14,4],[12,85,11],[5,15,91]])

sum = conf_arr.sum()

conf_arr = conf_arr * 100.0 / ( 1.0 * sum )

df_cm = pd.DataFrame(conf_arr, 
  index = [ 'Dog', 'Cat', 'Rabbit'],
  columns = ['Dog', 'Cat', 'Rabbit'])

fig = plt.figure()

plt.clf()

ax = fig.add_subplot(111)
ax.set_aspect(1)

cmap = sns.cubehelix_palette(light=1, as_cmap=True)

res = sn.heatmap(df_cm, annot=True, vmin=0.0, vmax=100.0, fmt='.2f', cmap=cmap)

res.invert_yaxis()

plt.yticks([0.5,1.5,2.5], [ 'Dog', 'Cat', 'Rabbit'],va='center')

plt.title('Confusion Matrix')

plt.savefig('confusion_matrix.png', dpi=100, bbox_inches='tight' )

plt.close()

Plot cldclass-lidar granule vertical profile using python 3

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC 
from pyhdf.HDF import *
from pyhdf.VS import *

import pprint
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import matplotlib.cm as cm

#----------------------------------------------------------------------------------------#

file_path = ''
file_name = '2008183012329_11573_CS_2B-CLDCLASS-LIDAR_GRANULE_P_R04_E02.hdf'

#----------------------------------------------------------------------------------------#
# Read HDF Files (VD data) Latitude & Longitude

f = HDF(file_path+file_name, SDC.READ) 
vs = f.vstart() 

Latitude = vs.attach('Latitude')
Longitude = vs.attach('Longitude')

a = Latitude[:]

Latitude.detach() 
Longitude.detach() 
vs.end()
f.close() 

#----------------------------------------------------------------------------------------#
# SDS Data

file = SD(file_path+file_name, SDC.READ)

file_info = file.info()
print(file_info)  # number of sds and metadata

print('---------- CloudLayerBase ----------')

sds_obj = file.select('CloudLayerBase') # select sds

CloudLayerBase = sds_obj.get() 

sds_info = sds_obj.info()

print(CloudLayerBase.shape)
print( sds_info )
print( sds_info[0], sds_info[1] )
print( 'sds attributes' )
pprint.pprint( sds_obj.attributes() )

for i in range(10):
	print([CloudLayerBase[i,j] for j in range(10)])

print('---------- CloudLayerTop ----------')

sds_obj = file.select('CloudLayerTop') # select sds

CloudLayerTop = sds_obj.get() 

sds_info = sds_obj.info()

print(CloudLayerTop.shape)
print( sds_info )
print( sds_info[0], sds_info[1] )
print( 'sds attributes' )
pprint.pprint( sds_obj.attributes() )

for i in range(10):
	print([CloudLayerTop[i,j] for j in range(10)])
	
print('---------- CloudPhase ----------')

sds_obj = file.select('CloudPhase') # select sds

CloudPhase = sds_obj.get() 

sds_info = sds_obj.info()

print(CloudPhase.shape)
print( sds_info )
print( sds_info[0], sds_info[1] )
print( 'sds attributes' )
pprint.pprint( sds_obj.attributes() )

for i in range(10):
	print([CloudPhase[i,j] for j in range(10)])


cldclass_lidar_start_idx = 000
cldclass_lidar_end_idx = 1000

tag = '2'
	
#----------------------------------------------------------------------------------------#
# Plot cldclass-lidar cloud layers

f = plt.figure()
ax = f.add_subplot(111)

#for i in range(CloudLayerBase.shape[0]):
for i in range(cldclass_lidar_start_idx,cldclass_lidar_end_idx):
	nb_cloud_layer = np.where(CloudLayerBase[i,:] < 0 )[0][0]
	for j in range(nb_cloud_layer):
		if CloudLayerBase[i,j] > 0 and CloudLayerTop[i,j] > 0.0:
			bar_xi = i
			bar_width = 1.0
			bar_base = CloudLayerBase[i,j] 
			bar_height = CloudLayerTop[i,j] - bar_base
			bar_color = '1.0'
			bar_edgecolor = '1.0'
			plt.bar(left=bar_xi, height=bar_height, width=bar_width, bottom=bar_base, color=bar_color, edgecolor=bar_edgecolor)

ax.set_facecolor('xkcd:lightblue')

plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='1.0', edgecolor='1.0', label='cloud')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='lightblue', edgecolor='lightblue', label='clear')

plt.legend()

plt.xlim(cldclass_lidar_start_idx,cldclass_lidar_end_idx)
plt.ylim(0,20)

plt.title('cldclass-lidar')

plt.savefig('plot_'+tag+'_cldclass_lidar_cloud_layers_method_01.png', dpi=100, bbox_inches='tight' )
#plt.show()
		
plt.close()

#----------------------------------------------------------------------------------------#
# Plot cldclass-lidar cloud layers + nb cloud layers (aggregation 1)

f = plt.figure()
ax = f.add_subplot(111)

#for i in range(CloudLayerBase.shape[0]):
for i in range(cldclass_lidar_start_idx,cldclass_lidar_end_idx):
	nb_cloud_layer = np.where(CloudLayerBase[i,:] < 0 )[0][0]
	for j in range(nb_cloud_layer):
		if CloudLayerBase[i,j] > 0 and CloudLayerTop[i,j] > 0.0:
			bar_xi = i
			bar_width = 1.0
			bar_base = CloudLayerBase[i,j] 
			bar_height = CloudLayerTop[i,j] - bar_base
			bar_color = '1.0'
			bar_edgecolor = '1.0'
			if nb_cloud_layer == 2:
				bar_color = 'blue'
				bar_edgecolor = 'blue'
			if nb_cloud_layer > 2:
				bar_color = 'salmon'
				bar_edgecolor = 'salmon'
			plt.bar(left=bar_xi, height=bar_height, width=bar_width, bottom=bar_base, color=bar_color, edgecolor=bar_edgecolor)

ax.set_facecolor('xkcd:lightblue')

plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='1.0', edgecolor='1.0', label='nb cloud layer = 1')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='blue', edgecolor='blue', label='nb cloud layer = 2')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='salmon', edgecolor='salmon', label='nb cloud layer > 2')

plt.legend()

plt.xlim(cldclass_lidar_start_idx,cldclass_lidar_end_idx)
plt.ylim(0,20)

plt.title('cldclass-lidar')

plt.savefig('plot_'+tag+'_cldclass_lidar_cloud_layers_and_nb_cloud_layers_method_01.png', dpi=100, bbox_inches='tight' )
#plt.show()
		
plt.close()

#----------------------------------------------------------------------------------------#
# Plot cldclass-lidar cloud layers + nb cloud layers (aggregation 2: merge two layers if
# same cloud phase and separation distance < threshold (3 for instance) )

f = plt.figure()
ax = f.add_subplot(111)

#for i in range(CloudLayerBase.shape[0]):
for i in range(cldclass_lidar_start_idx,cldclass_lidar_end_idx):
	nb_cloud_layer = np.where(CloudLayerBase[i,:] < 0 )[0][0]
	#for j in range(nb_cloud_layer):
	#	print(CloudLayerBase[i,j],CloudLayerTop[i,j])
	#print('-------')
	nb_cloud_layer_after_merging = nb_cloud_layer
	for j in range(nb_cloud_layer-1):	
			if ( CloudLayerBase[i,j+1] - CloudLayerTop[i,j] ) < 3.0:
				if CloudPhase[i,j+1] == CloudPhase[i,j]:
					nb_cloud_layer_after_merging = nb_cloud_layer_after_merging - 1
	for j in range(nb_cloud_layer):
		if CloudLayerBase[i,j] > 0 and CloudLayerTop[i,j] > 0.0:
			bar_xi = i
			bar_width = 1.0
			bar_base = CloudLayerBase[i,j] 
			bar_height = CloudLayerTop[i,j] - bar_base
			bar_color = '1.0'
			bar_edgecolor = '1.0'
			if nb_cloud_layer_after_merging == 2:
				bar_color = 'blue'
				bar_edgecolor = 'blue'
			if nb_cloud_layer_after_merging > 2:
				bar_color = 'salmon'
				bar_edgecolor = 'salmon'
			plt.bar(left=bar_xi, height=bar_height, width=bar_width, bottom=bar_base, color=bar_color, edgecolor=bar_edgecolor)

ax.set_facecolor('xkcd:lightblue')

plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='1.0', edgecolor='1.0', label='nb cloud layer = 1')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='blue', edgecolor='blue', label='nb cloud layer = 2')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='salmon', edgecolor='salmon', label='nb cloud layer > 2')

plt.legend()

plt.xlim(cldclass_lidar_start_idx,cldclass_lidar_end_idx)
plt.ylim(0,20)

plt.title('cldclass-lidar')

plt.savefig('plot_'+tag+'_cldclass_lidar_cloud_layers_and_nb_cloud_layers_aggregated.png', dpi=100, bbox_inches='tight' )
#plt.show()
		
plt.close()

#----------------------------------------------------------------------------------------#
# Plot cldclass-lidar cloud layers + nb cloud phase

f = plt.figure()
ax = f.add_subplot(111)

#for i in range(CloudLayerBase.shape[0]):
for i in range(cldclass_lidar_start_idx,cldclass_lidar_end_idx):
	nb_cloud_layer = np.where(CloudLayerBase[i,:] < 0 )[0][0]
	for j in range(nb_cloud_layer):
		if CloudLayerBase[i,j] > 0 and CloudLayerTop[i,j] > 0.0:
			bar_xi = i
			bar_width = 1.0
			bar_base = CloudLayerBase[i,j] 
			bar_height = CloudLayerTop[i,j] - bar_base
			bar_color = '1.0'
			bar_edgecolor = '1.0'
			if CloudPhase[i,j] == 3:
				bar_color = 'blue'
				bar_edgecolor = 'blue'
			if CloudPhase[i,j] == 2:
				bar_color = 'salmon'
				bar_edgecolor = 'salmon'			
			plt.bar(left=bar_xi, height=bar_height, width=bar_width, bottom=bar_base, color=bar_color, edgecolor=bar_edgecolor)

ax.set_facecolor('xkcd:lightblue')

plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='1.0', edgecolor='1.0', label='ice clouds')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='salmon', edgecolor='salmon', label='mixed clouds')
plt.bar(left=cldclass_lidar_start_idx, height=0, width=0, bottom=0, color='blue', edgecolor='blue', label='liquid clouds')

plt.legend()

plt.xlim(cldclass_lidar_start_idx,cldclass_lidar_end_idx)
plt.ylim(0,20)

plt.title('cldclass-lidar')

plt.savefig( 'plot_'+tag+'_cldclass_lidar_cloud_layers_and_thermodynamic_phase_method_01.png', dpi=100, bbox_inches='tight' )
#plt.show()
		
plt.close()

Python script modis cldclass-lidar analysis (main)

					
#!/usr/bin/env python

import numpy as np

import modis_cldclass_lidar_analysis_001 as analysis_001

#----------------------------------------------------------------------------------------#
# Settings

dict_settings = {}

dict_settings['input_path'] = ''

dict_settings['output_path'] = './media/images/modis_cldclass_lidar_analysis_main_py/'
dict_settings['output_file_name_prefix'] = 'modis_cldclass_lidar_analyisis_'

dict_settings['input_file_name'] = '2008_01_sample_01_modis_caliop_cldclass_lidar.txt'

data = np.loadtxt(
dict_settings['input_path'] + dict_settings['input_file_name'], skiprows=1)

fig_count = 0

#----------------------------------------------------------------------------------------#
# confusion matrices 

dict_settings['fig_title'] = 'modis c6 vs cldclass-lidar \n score = '
dict_settings['fig_ylabel'] = 'cldclass-lidar'
dict_settings['modis_cloud_optical_thickness_threshold'] = 4.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis-heidinger_test'] = False
dict_settings['product_source'] = 3

fig_count = analysis_001.multilayer_cloud_confusion_matrices(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 (without PH test) vs cldclass-lidar \n score = '
dict_settings['fig_ylabel'] = 'caliop 1km'
dict_settings['modis_cloud_optical_thickness_threshold'] = 4.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis_heidinger_test'] = True
dict_settings['product_source'] = 3

fig_count = analysis_001.multilayer_cloud_confusion_matrices(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 vs caliop 5km \n score = '
dict_settings['fig_ylabel'] = 'caliop 5km'
dict_settings['modis_cloud_optical_thickness_threshold'] = 4.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis-heidinger_test'] = False
dict_settings['product_source'] = 2

fig_count = analysis_001.multilayer_cloud_confusion_matrices(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'modis c6 vs caliop 1km \n score = '
dict_settings['fig_ylabel'] = 'caliop 1km'
dict_settings['modis_cloud_optical_thickness_threshold'] = 4.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis-heidinger_test'] = False
dict_settings['product_source'] = 1

fig_count = analysis_001.multilayer_cloud_confusion_matrices(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# confusion matrices (cloud phase)

dict_settings['fig_title'] = 'Cloud Thermodynamic Phase Comparisons \n modis c6 vs cldclass-lidar \n score = '
dict_settings['fig_ylabel'] = 'cldclass-lidar'
dict_settings['modis_cloud_optical_thickness_threshold'] = 0.0
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['modis_without_pavolonis-heidinger_test'] = False
dict_settings['product_source'] = 3

fig_count = analysis_001.cloud_thermodynamic_phase_confusion_matrices(data,dict_settings,fig_count)

#----------------------------------------------------------------------------------------#
# Cloud Layer Vertical Profile Histogram

dict_settings['fig_title'] = 'cldclass-lidar cloud layer vertical distribution'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['cloud_layer_logical_expression'] = 'int(data[i,18]) > 0'

fig_count = analysis_001.cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'cldclass-lidar cloud layer vertical distribution \n (nb cloud layer = 1)'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['cloud_layer_logical_expression'] = 'int(data[i,18]) == 1'

fig_count = analysis_001.cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'cldclass-lidar cloud layer vertical distribution \n (nb cloud layer = 2)'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['cloud_layer_logical_expression'] = 'int(data[i,18]) == 2'

fig_count = analysis_001.cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

dict_settings['fig_title'] = 'cldclass-lidar cloud layer vertical distribution \n (nb cloud layer > 1)'
dict_settings['modis_latitude_flag'] = True
dict_settings['modis_latitude_threshold'] = 60
dict_settings['cloud_layer_logical_expression'] = 'int(data[i,18]) > 1'

fig_count = analysis_001.cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count)

Python script modis cldclass-lidar analysis 001

					
#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sn
import pandas as pd

import seaborn as sns

#----------------------------------------------------------------------------------------#
# Confusion Matrices (Multilayer Clouds)

def multilayer_cloud_confusion_matrices(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]

	conf_arr = np.zeros((2,2))

	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]

		modis_cloud_optical_tickness = data[i,13] * 0.01
		
		if modis_cloud_optical_tickness > dict_settings['modis_cloud_optical_thickness_threshold']:
		
			if 'modis_without_pavolonis_heidinger_test' in dict_settings:
				if dict_settings['modis_without_pavolonis_heidinger_test']:
					if data[i,5] == 1: modis_multilayer_sds = modis_multilayer_sds - 3
		
			nb_cloud_layer = int( data[i,15] )			

			cldclass_lidar_phase_list = []
			cloud_layer_idx = 19
			for n in range(nb_cloud_layer):
				if nb_cloud_layer <= 10:
					if int(data[i,cloud_layer_idx+4]) == dict_settings['product_source']: # product source
						cldclass_lidar_phase_list.append(int(data[i,cloud_layer_idx+3]))
					cloud_layer_idx = cloud_layer_idx + 5			
			cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_phase_list)
		
			modis_multilayer_flag = -1
			if modis_multilayer_sds == 1 : modis_multilayer_flag = 0 
			if modis_multilayer_sds > 1: modis_multilayer_flag = 1
		
			cldclass_lidar_multilayer_flag = -1
			if cldclass_lidar_nb_cloud_layer == 1: cldclass_lidar_multilayer_flag = 0
			if cldclass_lidar_nb_cloud_layer > 1: cldclass_lidar_multilayer_flag = 1		
		
			if 	cldclass_lidar_multilayer_flag >= 0 and modis_multilayer_flag >= 0:
				conf_arr[modis_multilayer_flag,cldclass_lidar_multilayer_flag] = conf_arr[modis_multilayer_flag,cldclass_lidar_multilayer_flag] + 1
		
	sum = conf_arr.sum()
	
	score = ( conf_arr[0,0] + conf_arr[1,1] ) * 100.0 / sum
	conf_arr = conf_arr * 100.0 / ( 1.0 * sum )

	df_cm = pd.DataFrame(conf_arr, 
					  index = ['modis C6 monolayer cloud', 'modis C6 multilayer cloud'],
					  columns = [dict_settings['fig_ylabel']+'\n nb cloud layer = 1', dict_settings['fig_ylabel']+'\n nb cloud layer > 1'])

	fig = plt.figure()
	plt.clf()
	ax = fig.add_subplot(111)
	ax.set_aspect(1)

	cmap = sns.cubehelix_palette(light=1, as_cmap=True)

	res = sn.heatmap(df_cm, annot=True, vmin=0.0, vmax=100.0, fmt='.2f', cmap=cmap)
	
	res.invert_yaxis()
	
	plt.yticks([0.5,1.5], ['modis C6 \n multilayer cloud', 'modis C6 \n monolayer cloud'],va='center')

	plt.title(dict_settings['fig_title'] + str(round(score,2)))

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# Confusion Matrices (Cloud Thermodynamic Phase; cldclass-lidar monolayer)

def cloud_thermodynamic_phase_confusion_matrices(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]

	conf_arr = np.zeros((3,3))

	for i in range(data.shape[0]):
		
		modis_multilayer_sds = data[i,0]
		modis_cloud_phase_sds = data[i,6]		
		modis_cloud_optical_tickness = data[i,13] * 0.01
		
		if modis_cloud_optical_tickness > dict_settings['modis_cloud_optical_thickness_threshold']:
			
			nb_cloud_layer = int( data[i,15] )			

			cldclass_lidar_phase_list = []
			cloud_layer_idx = 19
			for n in range(nb_cloud_layer):
				if nb_cloud_layer <= 10:
					#print( data[i,cloud_layer_idx+4], data[i,cloud_layer_idx], data[i,cloud_layer_idx+1] )
					if int(data[i,cloud_layer_idx+4]) == 3: # source cldclass-lidar products
						cldclass_lidar_phase_list.append(int(data[i,cloud_layer_idx+3]))
					cloud_layer_idx = cloud_layer_idx + 5
			
			cldclass_lidar_nb_cloud_layer = len(cldclass_lidar_phase_list)
			
			#print(cldclass_lidar_phase_list,cldclass_lidar_nb_cloud_layer)
			#print('-------')
			
			#cldclass_lidar_multiple_phase = False
			#cloud_phase_layer_idx = 22
			#for j in range(int(cldclass_lidar_nb_cloud_layer)):
			#	if int(data[i,cloud_phase_layer_idx]) != int(cldclass_lidar_phase_sds):
			#		cldclass_lidar_multiple_phase = True
			#	cloud_phase_layer_idx = cloud_phase_layer_idx + 5
						
			if cldclass_lidar_nb_cloud_layer == 1:
			#if not cldclass_lidar_multiple_phase:
				
				cldclass_lidar_phase_sds = cldclass_lidar_phase_list[0]
				
				modis_cloud_phase_flag = -1
				if modis_cloud_phase_sds == 3: modis_cloud_phase_flag = 0 # ice
				if modis_cloud_phase_sds == 2: modis_cloud_phase_flag = 1 # liquid
				if modis_cloud_phase_sds == 4: modis_cloud_phase_flag = 2 # undetermined
		
				cldclass_lidar_cloud_phase_flag = -1
				if cldclass_lidar_phase_sds == 1: cldclass_lidar_cloud_phase_flag = 0 # ice
				if cldclass_lidar_phase_sds == 3: cldclass_lidar_cloud_phase_flag = 1 # liquid
				if cldclass_lidar_phase_sds == 2: cldclass_lidar_cloud_phase_flag = 2 # mixed	
		
				if 	modis_cloud_phase_flag >= 0 and cldclass_lidar_cloud_phase_flag >= 0:
					conf_arr[modis_cloud_phase_flag,cldclass_lidar_cloud_phase_flag] = conf_arr[modis_cloud_phase_flag,cldclass_lidar_cloud_phase_flag] + 1
		
	sum = conf_arr.sum()
	
	score = ( conf_arr[0,0] + conf_arr[1,1] ) * 100.0 / (conf_arr[0,0] + conf_arr[1,0] + conf_arr[0,1] + conf_arr[1,1])
		
	conf_arr = conf_arr * 100.0 / ( 1.0 * sum )

	df_cm = pd.DataFrame(conf_arr, 
					  index = ['modis C6 ice phase', 'modis C6 liquid phase', 'modis C6 undet. phase'],
					  columns = [dict_settings['fig_ylabel']+'\n ice phase', dict_settings['fig_ylabel']+'\n liquid phase', dict_settings['fig_ylabel']+'\n mixed phase'])

	fig = plt.figure()
	plt.clf()
	ax = fig.add_subplot(111)
	ax.set_aspect(1)

	cmap = sns.cubehelix_palette(light=1, as_cmap=True)

	res = sn.heatmap(df_cm, annot=True, vmin=0.0, vmax=100.0, fmt='.2f', cmap=cmap)
	
	res.invert_yaxis()
	
	plt.yticks([0.5,1.5,2.5], ['modis C6 \n undet. phase', 'modis C6 \n liquid phase', 'modis C6 \n ice phase'],va='center')

	plt.title(dict_settings['fig_title'] + str(round(score,2)))

	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

#----------------------------------------------------------------------------------------#
# Cloud Layer Vertical Distribution Histogram

np.set_printoptions(precision=2)

def cloud_layer_vertical_distribution_bar(data,dict_settings,fig_count):

	data = data[~np.all(data == 0, axis=1)]

	if 'modis_latitude_flag' in dict_settings:
		if dict_settings['modis_latitude_flag']:
			if 'modis_latitude_threshold' in dict_settings:				
				data = data[abs(data[:,7]) < dict_settings['modis_latitude_threshold']]
	
	cloud_layer_count_bar_x = [i*0.1 + 0.05 for i in range(200)]
	cloud_layer_count_bar_height = [0 for i in range(200)]
		
	for i in range(data.shape[0]):
	#for i in range(100):
		
		nb_cloud_layer = int( data[i,15] )
		
		cloud_layer_parameters = np.zeros((10,5))
		
		m = 19
		if nb_cloud_layer <= 10:
			for j in range(nb_cloud_layer):
				cloud_layer_parameters[j,0] = data[i,m]
				cloud_layer_parameters[j,1] = data[i,m+1]
				cloud_layer_parameters[j,2] = data[i,m+2]
				cloud_layer_parameters[j,3] = data[i,m+3]
				cloud_layer_parameters[j,4] = data[i,m+4]
				m = m + 5
			
			#print('nb_cloud_layer --->', nb_cloud_layer)
			#for j in range(nb_cloud_layer):
			#	print( [ "{0:0.1f}".format(value) for value in list(cloud_layer_parameters[j,:]) ] )
			#print('------')
		
			cloud_layer_caliop_count_bar = [0 for i in range(200)]

			'''
			for j in range(nb_cloud_layer):
				if cloud_layer_parameters[j,0] > 0.0 and cloud_layer_parameters[j,1] > 0.0:
					cloud_layer_top_idx = int(cloud_layer_parameters[j,0]/0.1)
					cloud_layer_base_idx = int(cloud_layer_parameters[j,1]/0.1)
					if int(cloud_layer_parameters[j,4]) == 3: # cldclass-lidar products
						for m in range(cloud_layer_base_idx,cloud_layer_top_idx+1):
							cloud_layer_count_bar_height[m] = cloud_layer_count_bar_height[m] + 1
			'''

			for j in range(nb_cloud_layer):
				if cloud_layer_parameters[j,0] > 0.0 and cloud_layer_parameters[j,1] > 0.0:
					cloud_layer_top_idx = int(cloud_layer_parameters[j,0]/0.1)
					cloud_layer_base_idx = int(cloud_layer_parameters[j,1]/0.1)
					if int(cloud_layer_parameters[j,4]) == 3: # cldclass-lidar products
						if eval(dict_settings['cloud_layer_logical_expression'] ):
							for m in range(cloud_layer_base_idx,cloud_layer_top_idx+1):
								cloud_layer_count_bar_height[m] = cloud_layer_count_bar_height[m] + 1
				
	
	plt.barh(cloud_layer_count_bar_x, cloud_layer_count_bar_height)
	
	plt.ylim(0,20)
	plt.xlim(0,50000)
	
	plt.title(dict_settings['fig_title'])
	
	plt.ylabel('altitude (km)')
	plt.xlabel('cloud layer counts')
		
	fig_count = fig_count + 1
	filename = dict_settings['output_file_name_prefix'] + '{:02d}'.format(fig_count) + '.png'
	plt.savefig(dict_settings['output_path'] + filename, dpi=100, bbox_inches='tight' )
	#plt.show()
	plt.close()

	return fig_count

About


To add a source code to your notebook just click on the link below. Sharing scientific source code allows to understand and reproduce your results:

Add a code

Search


Source codes can be retrieved using the following form: