Sample homogenization for a year of colocated MODIS and 2B-CLDCLASS-lidar data using python

					
#!/usr/bin/env python

from scipy.stats.kde import gaussian_kde

from random import uniform
from random import randint

import matplotlib.pyplot as plt
import numpy as np

import matplotlib.patches as mpatches

from statistics import mean

fig_count = 0

fig_ylim_max = 20000

path_to_media = "/Users/bmarcha1/Desktop/media/files/"

#----------------------------------------------------------------------------------------#
# Sample 1: concatenate monthly sample

data_1 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_01_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)
data_2 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_02_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)
data_3 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_03_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)
data_4 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_04_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)
data_5 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_05_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)
data_6 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_06_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)
data_7 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_07_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)
data_8 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_08_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)
data_9 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_09_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)
data_10 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_10_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)
data_11 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_11_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)
data_12 = np.loadtxt(path_to_media+'monthly_modis_cldclass_lidar/2008_12_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)

data = np.concatenate((data_1,data_2,data_3,data_4,data_5,data_6,
data_7,data_8,data_9,data_10,data_11,data_12), axis=0)

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

#----------------------------------------------------------------------------------------#
# Sample 1: plot zonal observations

nb_tot_obs = len(data[:,0])

bins = np.linspace(-90, 90, 180)

hist_counts_list, hist_x_bins_list, patches = plt.hist(data[:,0], bins=bins, color="blue", alpha=0.4)

plt.title('')
plt.xlabel('Latitudes')
plt.ylabel('Counts')

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

fig_count = fig_count + 1
filename ='2008_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
#plt.show()
plt.close()

#----------------------------------------------------------------------------------------#
# Sample 1: homogeneization (step 1)

data_obs = data[~np.all(data == 0, axis=1)]
data_cloudy = data_obs[data_obs[:,3] < 2]

nb_tot_obs = len(data_obs[:,0])
nb_tot_cloudy = len(data_cloudy[:,0])

print('nb_tot_obs: ',nb_tot_obs)
print('nb_tot_cloudy: ',nb_tot_cloudy)

bins = np.linspace(-90, 90, 180)

n_obs, bins_obs, patches_obs = plt.hist(data[:,0], bins=bins, color="#89bedc", alpha=1.0)
n_cloudy, bins_cloudy, patches_cloudy = plt.hist(data_cloudy[:,0], bins=bins, color="#0b559f", alpha=1.0)

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

plt.title('(a) before sample homogenization',fontsize=10)

plt.xlabel('Latitudes', fontsize=8)
plt.ylabel('Counts', fontsize=8)

plt.xticks(fontsize=8)
plt.yticks(fontsize=8)

pop_a = mpatches.Patch(color='#89bedc', label='Nb. of Observations')
pop_b = mpatches.Patch(color='#0b559f', label='Nb. of cloudy pixels')

plt.legend(handles=[pop_a,pop_b])

fig_count = fig_count + 1
filename ='2008_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
#plt.show()
plt.close()

prob_obs_bin = []
prob_cloudy_given_bin = []

for i in range(len(n_obs)):
	if n_obs[i] > 0.0:
		prob_obs_bin.append(1.0*n_obs[i]/nb_tot_obs)
		prob_cloudy_given_bin.append(n_cloudy[i]/n_obs[i])
	else:
		prob_obs_bin.append(0)
		prob_cloudy_given_bin.append(0)

prob_cloudy_bin = [prob_cloudy_given_bin[i]*prob_obs_bin[i]*nb_tot_obs/nb_tot_cloudy for i in range(len(prob_obs_bin))]

#----------------------------------------------------------------------------------------#
# Sample 1: homogeneization (step 2)

plt.bar([bins[i] for i in range(len(bins)-1)], [i*nb_tot_obs for i in prob_obs_bin])

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

plt.title('Random Sample Day Obs.') 
plt.xlabel('Latitude') 
plt.ylabel('Counts') 

fig_count = fig_count + 1
filename ='2008_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
plt.close()

homogenization_factor = 10.0 # ]0,100[
l = sorted([i for i in prob_obs_bin if i > 0.0])

homogenization_threshold_idx= int( len(l) * homogenization_factor / 100.0 )
homogenization_threshold = l[homogenization_threshold_idx]

for idx,i in enumerate(prob_obs_bin):
	if prob_obs_bin[idx] > homogenization_threshold:
		prob_obs_bin[idx] = homogenization_threshold

plt.bar([bins[i] for i in range(len(bins)-1)], [i*nb_tot_obs for i in prob_obs_bin])

plt.title('Random Sample Homogeneization') 
plt.xlabel('Latitude') 
plt.ylabel('Counts') 

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

fig_count = fig_count + 1
filename ='2008_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
plt.close()

plt.bar([bins[i] for i in range(len(bins)-1)], [i*nb_tot_cloudy for i in prob_cloudy_bin])

plt.title('Random Sample Day Cloudy Pixels') 
plt.xlabel('Latitude') 
plt.ylabel('Counts') 

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

fig_count = fig_count + 1
filename ='2008_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
plt.close()

prob_cloudy_bin = [prob_cloudy_given_bin[i]*prob_obs_bin[i]*nb_tot_obs/nb_tot_cloudy for i in range(len(prob_obs_bin))]

plt.bar([bins[i] for i in range(len(bins)-1)], [i*nb_tot_obs for i in prob_obs_bin], width=1.01, color="#89bedc")
plt.bar([bins[i] for i in range(len(bins)-1)], [i*nb_tot_cloudy for i in prob_cloudy_bin], width=1.01, color="#0b559f")

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

plt.title('(b) after sample homogenization',fontsize=10)

plt.xlabel('Latitudes', fontsize=8)
plt.ylabel('Counts', fontsize=8)

plt.xticks(fontsize=8)
plt.yticks(fontsize=8)

pop_a = mpatches.Patch(color='#89bedc', label='Nb. of Observations')
pop_b = mpatches.Patch(color='#0b559f', label='Nb. of cloudy pixels')

plt.legend(handles=[pop_a,pop_b], loc='upper center')



fig_count = fig_count + 1
filename ='2008_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
plt.close()

Sample homogenization for a month of colocated MODIS and 2B-CLDCLASS-lidar data using python

					
#!/usr/bin/env python

from scipy.stats.kde import gaussian_kde

from random import uniform
from random import randint

import matplotlib.pyplot as plt
import numpy as np

import matplotlib.patches as mpatches

from statistics import mean

fig_count = 0

fig_ylim_max = 2000

path_to_media = "/Users/bmarcha1/Desktop/media/files/"

#----------------------------------------------------------------------------------------#
# Sample 1

input_filename = path_to_media + 'monthly_modis_cldclass_lidar/2008_01_sample_00_modis_caliop_cldclass_lidar.txt'

data = np.loadtxt(input_filename)

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

nb_tot_obs = len(data[:,0])

#----------------------------------------------------------------------------------------#
# Sample 1: plot zonal observations

bins = np.linspace(-90, 90, 180)
	
hist_counts_list, hist_x_bins_list, patches = plt.hist(data[:,0], bins=bins, color="blue", alpha=0.4)

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

plt.title('zonal observations',fontsize=10)

plt.xlabel('Latitudes',fontsize=8)
plt.ylabel('Counts',fontsize=8)

fig_count = fig_count + 1
filename ='2008_01_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
#plt.show()
plt.close()

#----------------------------------------------------------------------------------------#
# Sample 1: homogeneization (step 1)

data_obs = data[~np.all(data == 0, axis=1)]
data_cloudy = data_obs[data_obs[:,3] < 2]

nb_tot_obs = len(data_obs[:,0])
nb_tot_cloudy = len(data_cloudy[:,0])

print('nb_tot_obs: ',nb_tot_obs)
print('nb_tot_cloudy: ',nb_tot_cloudy)

bins = np.linspace(-90, 90, 180)
	
n_obs, bins_obs, patches_obs = plt.hist(data[:,0], bins=bins, color="#89bedc", alpha=1.0)
n_cloudy, bins_cloudy, patches_cloudy = plt.hist(data_cloudy[:,0], bins=bins, color="#0b559f", alpha=1.0)

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

plt.title('Jamuary 2008 colocated MODIS-2B-CLDCLASS-lidar sample \n (a) before sample homogenization',fontsize=10)

plt.xlabel('Latitudes', fontsize=8)
plt.ylabel('Counts', fontsize=8)

plt.xticks(fontsize=8)
plt.yticks(fontsize=8)

pop_a = mpatches.Patch(color='#89bedc', label='Nb. of Observations')
pop_b = mpatches.Patch(color='#0b559f', label='Nb. of cloudy pixels')

plt.legend(handles=[pop_a,pop_b])

fig_count = fig_count + 1
filename ='2008_01_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
#plt.show()
plt.close()

prob_obs_bin = []
prob_cloudy_given_bin = []

for i in range(len(n_obs)):
	if n_obs[i] > 0.0:
		prob_obs_bin.append(1.0*n_obs[i]/nb_tot_obs)
		prob_cloudy_given_bin.append(n_cloudy[i]/n_obs[i])
	else:
		prob_obs_bin.append(0)
		prob_cloudy_given_bin.append(0)
	
prob_cloudy_bin = [prob_cloudy_given_bin[i]*prob_obs_bin[i]*nb_tot_obs/nb_tot_cloudy for i in range(len(prob_obs_bin))]

#----------------------------------------------------------------------------------------#
# Sample 1: homogeneization (step 2)

plt.bar([bins[i] for i in range(len(bins)-1)], [i*nb_tot_obs for i in prob_obs_bin], color="royalblue")

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

fig_count = fig_count + 1
filename ='2008_01_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
plt.close()

homogenization_factor = 30.0 # ]0,100[
l = sorted([i for i in prob_obs_bin if i > 0.0])

homogenization_threshold_idx= int( len(l) * homogenization_factor / 100.0 )
homogenization_threshold = l[homogenization_threshold_idx]

for idx,i in enumerate(prob_obs_bin):
	if prob_obs_bin[idx] > homogenization_threshold:
		prob_obs_bin[idx] = homogenization_threshold

plt.bar([bins[i] for i in range(len(bins)-1)], [i*nb_tot_obs for i in prob_obs_bin], color="royalblue")

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

fig_count = fig_count + 1
filename ='2008_01_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
plt.close()

plt.bar([bins[i] for i in range(len(bins)-1)], [i*nb_tot_cloudy for i in prob_cloudy_bin], color="lightcoral")

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

fig_count = fig_count + 1
filename ='2008_01_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
plt.close()

prob_cloudy_bin = [prob_cloudy_given_bin[i]*prob_obs_bin[i]*nb_tot_obs/nb_tot_cloudy for i in range(len(prob_obs_bin))]

plt.bar([bins[i] for i in range(len(bins)-1)], [i*nb_tot_obs for i in prob_obs_bin], width=1.01, color="#89bedc")
plt.bar([bins[i] for i in range(len(bins)-1)], [i*nb_tot_cloudy for i in prob_cloudy_bin], width=1.01, color="#0b559f")

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

plt.title('Jamuary 2008 colocated MODIS-2B-CLDCLASS-lidar sample \n (b) after sample homogenization',fontsize=10)

plt.xlabel('Latitudes', fontsize=8)
plt.ylabel('Counts', fontsize=8)

plt.xticks(fontsize=8)
plt.yticks(fontsize=8)

pop_a = mpatches.Patch(color='#89bedc', label='Nb. of Observations')
pop_b = mpatches.Patch(color='#0b559f', label='Nb. of cloudy pixels')

plt.legend(handles=[pop_a,pop_b], loc='upper center')

fig_count = fig_count + 1
filename ='2008_01_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
plt.close()

#----------------------------------------------------------------------------------------#
# Sample 2: homogenization

input_filename = path_to_media + 'monthly_modis_cldclass_lidar/2008_01_sample_01_modis_caliop_cldclass_lidar.txt'

sample_2_data = np.loadtxt(input_filename)
sample_2_data = sample_2_data[~np.all(sample_2_data == 0, axis=1)]

nb_tot_cloudy= len(sample_2_data[:,7])

bins = np.linspace(-90, 90, 180)
	
hist_counts_list, hist_x_bins_list, patches = plt.hist(sample_2_data[:,7], bins=bins, color="lightcoral", alpha=0.4)

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

plt.title('')

plt.xlabel('Latitudes')
plt.ylabel('Counts')

fig_count = fig_count + 1
filename ='2008_01_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
plt.close()

plt.bar([bins[i] for i in range(len(bins)-1)], [i for i in hist_counts_list], color="lightcoral")

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

fig_count = fig_count + 1
filename ='2008_01_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
plt.close()

n_bin_max_list = [i*nb_tot_cloudy for i in prob_cloudy_bin]

plt.bar([bins[i] for i in range(len(bins)-1)], [i if i < n_bin_max_list[idx] else n_bin_max_list[idx] for idx,i in enumerate(hist_counts_list)], color="lightcoral")

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

fig_count = fig_count + 1
filename ='2008_01_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
plt.close()

data_2 = np.zeros(sample_2_data.shape)
nb_points = 0
for i in range(len(hist_counts_list)):
	if hist_counts_list[i] > 0.0:
		data_bin = sample_2_data[(sample_2_data[:,7] >= bins[i]) & (sample_2_data[:,7] <  bins[i+1]) ]
		np.random.shuffle(data_bin)
		n_bin_count = 0
		for j in range(data_bin.shape[0]):
			if n_bin_count < n_bin_max_list[i]:
				data_2[nb_points,:] = data_bin[j,:]
				nb_points = nb_points + 1
				n_bin_count = n_bin_count + 1
data_2 = data_2[~np.all(data_2 == 0, axis=1)]

n, bins, patches = plt.hist(data_2[:,7], bins=bins, color="lightcoral", alpha=0.4)

plt.xlim(-90.0,90.0)
plt.ylim(0,fig_ylim_max)
plt.grid()

plt.title('Homogenized sample')

plt.xlabel('Latitudes')
plt.ylabel('Counts')

fig_count = fig_count + 1
filename ='2008_01_sample_homogenization_'+'{:02d}'.format(fig_count)+'.png'
plt.savefig(filename, dpi=100, bbox_inches='tight' )
#plt.show()
plt.close()

How to plot scatter points on a global map using matplotlib and basemap ?

					
#!/usr/bin/env python

from mpl_toolkits.basemap import Basemap, cm

from scipy import stats, mgrid, c_, reshape, random, rot90
from pylab import *
from mpl_toolkits.axes_grid1 import make_axes_locatable

import numpy as np
import seaborn as sn
import pandas as pd

import matplotlib as mpl

mpl.style.use('seaborn')

data = np.loadtxt('2008_01_sample_00_modis_caliop_cldclass_lidar.txt', skiprows=1)

np.random.shuffle(data)

data = data[:10000,:]

#---------------------------------------------------------------------------------------#
# Scatter Plots Observations

fig = plt.figure()

ax = fig.add_subplot(111)

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.))

m.scatter(data[:,1],data[:,0],s=2)

plt.title('Random sample with a reservoir for January 2008',fontsize='10')

plt.xlim(-180,180)
plt.ylim(-90,90)

plt.xlabel('Longitude',fontsize='8')
plt.ylabel('Latitude',fontsize='8')

filename ='random_sample_global_map.png'
plt.savefig(filename, dpi=200, bbox_inches='tight' )
plt.show()
#plt.close()

Plot MODIS C6 L2 Cloud Mask 1km with python and matplotlib (version 2)

					
#!/usr/bin/env python

#Author: Benjamin Marchant

from pyhdf.SD import SD, SDC 

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

#----------------------------------------------------------------------------------------#
# Inputs

file_name = 'MYD06_L2.A2008015.1435.006.2013342100940.hdf'

#----------------------------------------------------------------------------------------#
# Read HDF Files

file = SD(file_name, SDC.READ)

data_selected_id = file.select('Cloud_Mask_1km')

data = data_selected_id.get()

data_shape = data.shape

#----------------------------------------------------------------------------------------#
# Read bits function

def bits_stripping(bit_start,bit_count,value):
	bitmask=pow(2,bit_start+bit_count)-1
	return np.right_shift(np.bitwise_and(value,bitmask),bit_start)

#----------------------------------------------------------------------------------------#
# Plot data

#status_flag = bits_stripping(0,1,data[i,j,0]) 

cloud_mask_flag = bits_stripping(1,2,data[:,:,0])

print(cloud_mask_flag[100,:])

cmap =  [(1.0,1.0,1.0)] + [(1.0, 0.0, 0.0)] + [(65.0/255,105.0/255,225.0/255)] + [(0.0,0.0,0.0)]
cmap = mpl.colors.ListedColormap(cmap)

bounds = [-0.5, 0.5, 1.5, 2.5, 3.5]

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

img = plt.imshow(np.fliplr(cloud_mask_flag), cmap=cmap, norm=norm,
				  interpolation='none', origin='lower')

cbar_bounds = [-0.5, 0.5, 1.5, 2.5, 3.5]
cbar_ticks = [ 0, 1, 2, 3]  
cbar_labels = ['Confident Cloudy', 'Probably Cloudy','Probably Clear ','Confident Clear']  
		   
cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_labels, fontsize=8)

plt.title('Cloud Mask 1km \n MYD06 C6 (2008-01-15; 14h35)', fontsize=8)


l = [int(i) for i in np.linspace(0,data_shape[1],6)]
plt.xticks(l, [i for i in reversed(l)], rotation=0, fontsize=8 )

l = [int(i) for i in np.linspace(0,data_shape[0],9)]
plt.yticks(l, l, rotation=0, fontsize=8 )

plt.xticks(fontsize=8)
plt.yticks(fontsize=8)

plt.savefig("myd06_c6_cloud_mask_1km_baseline_v2.png", 
			bbox_inches='tight', dpi=100)
#plt.show()

Exemple de code python pour créer un diagramme de Voronoï

					
import matplotlib.pyplot as plt
import numpy as np 
import math

x_min, x_max = 0.0, 10.0
y_min, y_max = 0.0, 10.0

x_pts = [2.0,3.8,5.0,7.7,8.0]
y_pts = [2.0,8.2,5.0,2.4,7.6]

plt.scatter(x_pts,y_pts, color='gray')

x = np.linspace(0, 10, 100)

for idx in range(len(x_pts)):

	voronoi_segment_list = []

	#----------------------------------------------------------------------------------------#
	# Step 1: find nearest point

	if idx < len(x_pts)-1 :
		jdx_min = idx + 1
	else:
		jdx_min = idx - 1

	d_min = math.sqrt( (x_pts[idx]-x_pts[jdx_min])**2 + (y_pts[idx]-y_pts[jdx_min])**2 )

	for jdx,pts in enumerate(x_pts):
		if jdx != idx:
			d = math.sqrt( (x_pts[idx]-x_pts[jdx])**2 + (y_pts[idx]-y_pts[jdx])**2 )
			if d < d_min:
				d_min = d
				jdx_min = jdx
				
	#----------------------------------------------------------------------------------------#
	# Step 2: find all perpendicular lines

	perp_line_list = []

	for jdx,pts in enumerate(x_pts):
		if jdx != idx:
		
			x_middle = 0.5 * ( x_pts[idx] + x_pts[jdx] )
			y_middle = 0.5 * ( y_pts[idx] + y_pts[jdx] )
		
			slope = ( y_pts[jdx] - y_pts[idx] ) / ( x_pts[jdx] - x_pts[idx] )
		
			perp_slope = -1.0 / slope 
			perp_intercept = y_middle - perp_slope * x_middle

			if jdx == jdx_min: first_perp_line_idx = len(perp_line_list)	
		
			perp_line_list.append((perp_slope,perp_intercept,x_middle,y_middle))
		
			#y = perp_slope * x + perp_intercept

			#plt.plot(x,y, color='gray', linestyle='--')

	#----------------------------------------------------------------------------------------#
	# Step 3: find the first voronoi segment

	x_middle_0 = 0.5 * ( x_pts[idx] + x_pts[jdx_min] )
	y_middle_0 = 0.5 * ( y_pts[idx] + y_pts[jdx_min] )

	slope_0 = perp_line_list[first_perp_line_idx][0]
	intercept_0 = perp_line_list[first_perp_line_idx][1]

	#y = slope_0 * x + intercept_0

	#plt.plot(x,y, color='red', linestyle='--')

	intersection_list = []

	intersection_list.append((x_min,slope_0 * x_min + intercept_0,slope_0,intercept_0))

	intersection_list.append((x_max,slope_0 * x_max + intercept_0,slope_0,intercept_0))

	for line in perp_line_list:

		slope = line[0]
		intercept = line[1]

		if slope != slope_0 and intercept != intercept_0:
	
			intersection_x = (intercept - intercept_0) / ( slope_0 - slope )

			intersection_list.append((intersection_x,slope * intersection_x + intercept,slope,intercept))

	intersection_list.sort(key=lambda x: x[0])

	#print(intersection_list)
	#print(' ')

	line_0_list = []

	for i in range(len(intersection_list)-1):
	
		if x_middle_0>intersection_list[i][0] and x_middle_0<intersection_list[i+1][0]:
		
			voronoi_segment_list.append((
			intersection_list[i][0],
			intersection_list[i][1],
			intersection_list[i+1][0],
			intersection_list[i+1][1],
			slope_0,
			intercept_0))
		
			line_0_list.append((
			intersection_list[i][0],
			intersection_list[i][1],
			intersection_list[i][2],
			intersection_list[i][3],))

			line_0_list.append((
			intersection_list[i+1][0],
			intersection_list[i+1][1],
			intersection_list[i+1][2],
			intersection_list[i+1][3],))

	#print('voronoi_segment_list',voronoi_segment_list)
	#print(' ')

	#print('line_0_list',line_0_list)
	#print(' ')

	#----------------------------------------------------------------------------------------#
	# Step 4: find all voronoi segments

	while line_0_list:

		new_line_0_list = []
		new_voronoi_segment_list = []
	
		for line_0 in line_0_list:

			slope_0 = line_0[2]
			intercept_0 = line_0[3]
		
			x_0 = line_0[0]
			y_0 = line_0[1]
		
			#print('slope_0,intercept_0',slope_0,intercept_0,x_0,y_0)
			#print(' ')
		
			#y = slope_0 * x + intercept_0

			#plt.plot(x,y, color='yellow', linestyle='--')

			intersection_list = []
		
			fx_min = slope_0 * x_min + intercept_0
			d = math.sqrt( (x_min-x_0)**2 + (fx_min-y_0)**2)
			intersection_list.append((x_min,fx_min,d,slope_0,intercept_0))

			fx_min = slope_0 * x_max + intercept_0
			d = math.sqrt( (x_min-x_0)**2 + (fx_min-y_0)**2)
			intersection_list.append((x_max,fx_min,d,slope_0,intercept_0))

			for line in perp_line_list:

				slope = line[0]
				intercept = line[1]
			
				if slope != slope_0 and intercept != intercept_0:
	
					intersection_x = (intercept - intercept_0) / ( slope_0 - slope )
				
					fx_min = slope * intersection_x + intercept
					d = math.sqrt( (intersection_x-x_0)**2 + (fx_min-y_0)**2)
				
					intersection_list.append((intersection_x,fx_min,d,slope,intercept))

			intersection_list.sort(key=lambda x: x[0])

			#print('intersection_list',intersection_list)
			#print(' ')	
		
			intersection_list_filtered = []
		
			for i,intersection in enumerate(intersection_list):
		
				#if intersection[2] != 0.0: 
				if abs(intersection[2]) > 0.000000001: 
			
					for voronoi_segment in voronoi_segment_list:
				
						a = y_pts[idx] - ( voronoi_segment[4] * x_pts[idx]  + voronoi_segment[5] )
						b = intersection[1] - ( voronoi_segment[4] * intersection[0]  + voronoi_segment[5] )
						
						if a > 0.0 and b > 0.0: intersection_list_filtered.append(intersection)
						if a < 0.0 and b < 0.0: intersection_list_filtered.append(intersection)

		
			intersection_list_filtered.sort(key=lambda x: x[2])
		
			#print('intersection_list filtered',intersection_list_filtered)
			#print(' ')	
			#print('----------------')		
		
		
			intersection_slope = intersection_list_filtered[0][3]
			intersection_intercept = intersection_list_filtered[0][4]

			cond = True		
			for voronoi_segment in voronoi_segment_list:		
				if intersection_slope == voronoi_segment[4] and intersection_intercept == voronoi_segment[5]:
					cond = False
		
			for new_line_0 in new_line_0_list:		
				if intersection_slope == new_line_0[2] and intersection_intercept == new_line_0[3]:
					cond = False
		
			if cond:
				new_line_0_list.append((
				intersection_list_filtered[0][0],
				intersection_list_filtered[0][1],
				intersection_list_filtered[0][3],
				intersection_list_filtered[0][4]))
		
			new_voronoi_segment_list.append((
			x_0,
			y_0,
			intersection_list_filtered[0][0],
			intersection_list_filtered[0][1],
			slope_0,
			intercept_0))

		for new_voronoi_segment in new_voronoi_segment_list:
			voronoi_segment_list.append(new_voronoi_segment)
	
		#print('new_line_0_list',new_line_0_list)
		#print(' ')		

		line_0_list = new_line_0_list

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

	#print(voronoi_segment_list,len(voronoi_segment_list))

	for voronoi_segment in voronoi_segment_list:
	
		voronoi_segment_x = [voronoi_segment[0],voronoi_segment[2]]
		voronoi_segment_y = [voronoi_segment[1],voronoi_segment[3]]

		plt.plot(voronoi_segment_x, voronoi_segment_y, color='coral', linestyle='-')

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

plt.title("How to create a Voronoi diagram \n based on euclidean distance using python ?", fontsize=10)
plt.xlabel('x',fontsize=8)
plt.ylabel('y',fontsize=8)

plt.xlim(0,10)
plt.ylim(0,10)

plt.grid()

plt.savefig("voronoi_01.png")
#plt.show()

A python code to create Voronoi diagram based on Euclidean distance

					
import matplotlib.pyplot as plt
import numpy as np 
import math

x_min, x_max = 0.0, 10.0
y_min, y_max = 0.0, 10.0

x_pts = [2.0,3.8,5.0,7.7,8.0]
y_pts = [2.0,8.2,5.0,2.4,7.6]

plt.scatter(x_pts,y_pts, color='gray')

x = np.linspace(0, 10, 100)

for idx in range(len(x_pts)):

	voronoi_segment_list = []

	#----------------------------------------------------------------------------------------#
	# Step 1: find nearest point

	if idx < len(x_pts)-1 :
		jdx_min = idx + 1
	else:
		jdx_min = idx - 1

	d_min = math.sqrt( (x_pts[idx]-x_pts[jdx_min])**2 + (y_pts[idx]-y_pts[jdx_min])**2 )

	for jdx,pts in enumerate(x_pts):
		if jdx != idx:
			d = math.sqrt( (x_pts[idx]-x_pts[jdx])**2 + (y_pts[idx]-y_pts[jdx])**2 )
			if d < d_min:
				d_min = d
				jdx_min = jdx
				
	#----------------------------------------------------------------------------------------#
	# Step 2: find all perpendicular lines

	perp_line_list = []

	for jdx,pts in enumerate(x_pts):
		if jdx != idx:
		
			x_middle = 0.5 * ( x_pts[idx] + x_pts[jdx] )
			y_middle = 0.5 * ( y_pts[idx] + y_pts[jdx] )
		
			slope = ( y_pts[jdx] - y_pts[idx] ) / ( x_pts[jdx] - x_pts[idx] )
		
			perp_slope = -1.0 / slope 
			perp_intercept = y_middle - perp_slope * x_middle

			if jdx == jdx_min: first_perp_line_idx = len(perp_line_list)	
		
			perp_line_list.append((perp_slope,perp_intercept,x_middle,y_middle))
		
			#y = perp_slope * x + perp_intercept

			#plt.plot(x,y, color='gray', linestyle='--')

	#----------------------------------------------------------------------------------------#
	# Step 3: find the first voronoi segment

	x_middle_0 = 0.5 * ( x_pts[idx] + x_pts[jdx_min] )
	y_middle_0 = 0.5 * ( y_pts[idx] + y_pts[jdx_min] )

	slope_0 = perp_line_list[first_perp_line_idx][0]
	intercept_0 = perp_line_list[first_perp_line_idx][1]

	#y = slope_0 * x + intercept_0

	#plt.plot(x,y, color='red', linestyle='--')

	intersection_list = []

	intersection_list.append((x_min,slope_0 * x_min + intercept_0,slope_0,intercept_0))

	intersection_list.append((x_max,slope_0 * x_max + intercept_0,slope_0,intercept_0))

	for line in perp_line_list:

		slope = line[0]
		intercept = line[1]

		if slope != slope_0 and intercept != intercept_0:
	
			intersection_x = (intercept - intercept_0) / ( slope_0 - slope )

			intersection_list.append((intersection_x,slope * intersection_x + intercept,slope,intercept))

	intersection_list.sort(key=lambda x: x[0])

	#print(intersection_list)
	#print(' ')

	line_0_list = []

	for i in range(len(intersection_list)-1):
	
		if x_middle_0>intersection_list[i][0] and x_middle_0<intersection_list[i+1][0]:
		
			voronoi_segment_list.append((
			intersection_list[i][0],
			intersection_list[i][1],
			intersection_list[i+1][0],
			intersection_list[i+1][1],
			slope_0,
			intercept_0))
		
			line_0_list.append((
			intersection_list[i][0],
			intersection_list[i][1],
			intersection_list[i][2],
			intersection_list[i][3],))

			line_0_list.append((
			intersection_list[i+1][0],
			intersection_list[i+1][1],
			intersection_list[i+1][2],
			intersection_list[i+1][3],))

	#print('voronoi_segment_list',voronoi_segment_list)
	#print(' ')

	#print('line_0_list',line_0_list)
	#print(' ')

	#----------------------------------------------------------------------------------------#
	# Step 4: find all voronoi segments

	while line_0_list:

		new_line_0_list = []
		new_voronoi_segment_list = []
	
		for line_0 in line_0_list:

			slope_0 = line_0[2]
			intercept_0 = line_0[3]
		
			x_0 = line_0[0]
			y_0 = line_0[1]
		
			#print('slope_0,intercept_0',slope_0,intercept_0,x_0,y_0)
			#print(' ')
		
			#y = slope_0 * x + intercept_0

			#plt.plot(x,y, color='yellow', linestyle='--')

			intersection_list = []
		
			fx_min = slope_0 * x_min + intercept_0
			d = math.sqrt( (x_min-x_0)**2 + (fx_min-y_0)**2)
			intersection_list.append((x_min,fx_min,d,slope_0,intercept_0))

			fx_min = slope_0 * x_max + intercept_0
			d = math.sqrt( (x_min-x_0)**2 + (fx_min-y_0)**2)
			intersection_list.append((x_max,fx_min,d,slope_0,intercept_0))

			for line in perp_line_list:

				slope = line[0]
				intercept = line[1]
			
				if slope != slope_0 and intercept != intercept_0:
	
					intersection_x = (intercept - intercept_0) / ( slope_0 - slope )
				
					fx_min = slope * intersection_x + intercept
					d = math.sqrt( (intersection_x-x_0)**2 + (fx_min-y_0)**2)
				
					intersection_list.append((intersection_x,fx_min,d,slope,intercept))

			intersection_list.sort(key=lambda x: x[0])

			#print('intersection_list',intersection_list)
			#print(' ')	
		
			intersection_list_filtered = []
		
			for i,intersection in enumerate(intersection_list):
		
				#if intersection[2] != 0.0: 
				if abs(intersection[2]) > 0.000000001: 
			
					for voronoi_segment in voronoi_segment_list:
				
						a = y_pts[idx] - ( voronoi_segment[4] * x_pts[idx]  + voronoi_segment[5] )
						b = intersection[1] - ( voronoi_segment[4] * intersection[0]  + voronoi_segment[5] )
						
						if a > 0.0 and b > 0.0: intersection_list_filtered.append(intersection)
						if a < 0.0 and b < 0.0: intersection_list_filtered.append(intersection)

		
			intersection_list_filtered.sort(key=lambda x: x[2])
		
			#print('intersection_list filtered',intersection_list_filtered)
			#print(' ')	
			#print('----------------')		
		
		
			intersection_slope = intersection_list_filtered[0][3]
			intersection_intercept = intersection_list_filtered[0][4]

			cond = True		
			for voronoi_segment in voronoi_segment_list:		
				if intersection_slope == voronoi_segment[4] and intersection_intercept == voronoi_segment[5]:
					cond = False
		
			for new_line_0 in new_line_0_list:		
				if intersection_slope == new_line_0[2] and intersection_intercept == new_line_0[3]:
					cond = False
		
			if cond:
				new_line_0_list.append((
				intersection_list_filtered[0][0],
				intersection_list_filtered[0][1],
				intersection_list_filtered[0][3],
				intersection_list_filtered[0][4]))
		
			new_voronoi_segment_list.append((
			x_0,
			y_0,
			intersection_list_filtered[0][0],
			intersection_list_filtered[0][1],
			slope_0,
			intercept_0))

		for new_voronoi_segment in new_voronoi_segment_list:
			voronoi_segment_list.append(new_voronoi_segment)
	
		#print('new_line_0_list',new_line_0_list)
		#print(' ')		

		line_0_list = new_line_0_list

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

	#print(voronoi_segment_list,len(voronoi_segment_list))

	for voronoi_segment in voronoi_segment_list:
	
		voronoi_segment_x = [voronoi_segment[0],voronoi_segment[2]]
		voronoi_segment_y = [voronoi_segment[1],voronoi_segment[3]]

		plt.plot(voronoi_segment_x, voronoi_segment_y, color='coral', linestyle='-')

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

plt.title("How to create a Voronoi diagram \n based on euclidean distance using python ?", fontsize=10)
plt.xlabel('x',fontsize=8)
plt.ylabel('y',fontsize=8)

plt.xlim(0,10)
plt.ylim(0,10)

plt.grid()

plt.savefig("voronoi_01.png")
#plt.show()

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)

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