# How to plot MODIS cloud re and tau LUT (Nakajima and King's diagram) using python ?

Ben June 27, 2019

An example of how to plot MODIS cloud re and tau LUT (Nakajima and King's plot) in python using matplotlib:

### How to get and read the Data

Note: the LUTs are available here

from pyhdf.SD import SD, SDC
from scipy import interpolate

import matplotlib.pyplot as plt
import numpy as np
import math
import matplotlib.patches as mpatches
import matplotlib.cm as cm

print( file.info() )

#----------------------------------------------------------------------------------------#
# Get SDS info

datasets_dic = file.datasets()

for idx,sds in enumerate(datasets_dic.keys()):
print( idx,sds )
sds_obj = file.select(sds)
print( sds_obj.info() )
print()
print()

#----------------------------------------------------------------------------------------#
# Get the Data

reflectanceWater = file.select('reflectanceWater')
reflectanceIce = file.select('reflectanceIce')
WaveLengths = file.select('WaveLengths')
OpticalThickness = file.select('OpticalThickness')

reflectanceWater_data = reflectanceWater.get()
reflectanceIce_data = reflectanceIce.get()
WaveLengths_data = WaveLengths.get()
OpticalThickness_data = OpticalThickness.get()

liq_lut_shape = reflectanceWater_data.shape
ice_lut_shape = reflectanceIce_data.shape

### Nakajima and King's diagram

# band idx:
# 0 --> 0.65899998
# 1 --> 0.86500001
# 2 --> 1.24000001
# 3 --> 1.63999999
# 4 --> 2.13000011
# 5 --> 3.75

x_band_idx = 0
y_band_idx = 4

phase = 'ice' # liquid or ice

fig = plt.figure()

if phase == 'liquid':
x = [reflectanceWater_data[index_re,x_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]
y = [reflectanceWater_data[index_re,y_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]

if phase == 'ice':
x = [reflectanceIce_data[index_re,x_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]
y = [reflectanceIce_data[index_re,y_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]

plt.plot(x,y, 'steelblue',linewidth=0.5)

if phase == 'liquid':
x = [reflectanceWater_data[:,x_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]
y = [reflectanceWater_data[:,y_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]

if phase == 'ice':
x = [reflectanceIce_data[:,x_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]
y = [reflectanceIce_data[:,y_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]

plt.plot(x,y, 'coral',linewidth=0.5)

pop_a = mpatches.Patch(color='coral', label=r'Cloud Effective Radius ($\mu m$)')
pop_b = mpatches.Patch(color='steelblue', label='Cloud Optical Thickness')

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

if phase == 'liquid': plt.title("MODIS Liquid Cloud LUT",fontsize=8)
if phase == 'ice': plt.title("MODIS Ice Cloud LUT",fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[x_band_idx]) + r'\mu m$'
plt.xlabel(label,fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[y_band_idx]) + r'\mu m$'
plt.ylabel(label,fontsize=8)

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

plt.show()
plt.close()

### Add Re and Tau labels

#print(reflectanceIce_data)
#print(OpticalThickness_data)

x_band_idx = 0
y_band_idx = 4

phase = 'ice' # liquid or ice

fig = plt.figure()

if phase == 'liquid':
x = [reflectanceWater_data[index_re,x_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]
y = [reflectanceWater_data[index_re,y_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]

x_label_re = [reflectanceWater_data[index_re,x_band_idx,liq_lut_shape[2]-1] for index_re in np.arange(liq_lut_shape[0]) ]
y_label_re = [reflectanceWater_data[index_re,y_band_idx,liq_lut_shape[2]-1] for index_re in np.arange(liq_lut_shape[0]) ]

if phase == 'ice':
x = [reflectanceIce_data[index_re,x_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]
y = [reflectanceIce_data[index_re,y_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]

x_label_re = [reflectanceIce_data[index_re,x_band_idx,ice_lut_shape[2]-1] for index_re in np.arange(ice_lut_shape[0]) ]
y_label_re = [reflectanceIce_data[index_re,y_band_idx,ice_lut_shape[2]-1] for index_re in np.arange(ice_lut_shape[0]) ]

plt.plot(x,y, 'steelblue',linewidth=0.5)

for i,j,z in zip(x_label_re, y_label_re,label_re):
plt.text( i+0.04, j, z, ha='center', va='center', fontsize=8, color='coral')

if phase == 'liquid':
x = [reflectanceWater_data[:,x_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]
y = [reflectanceWater_data[:,y_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]

idx_tau_list = [0,3,6,9,12,15,18,21,34]

x_label = [reflectanceWater_data[0,x_band_idx,index_tau] for index_tau in idx_tau_list ]
y_label = [reflectanceWater_data[0,y_band_idx,index_tau] for index_tau in idx_tau_list ]

label = [str(OpticalThickness_data[re_idx]) for re_idx in idx_tau_list ]

if phase == 'ice':
x = [reflectanceIce_data[:,x_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]
y = [reflectanceIce_data[:,y_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]

idx_tau_list = [0,3,6,9,12,15,18,21,34]

x_label = [reflectanceIce_data[0,x_band_idx,index_tau] for index_tau in idx_tau_list ]
y_label = [reflectanceIce_data[0,y_band_idx,index_tau] for index_tau in idx_tau_list ]

label = [str(OpticalThickness_data[re_idx]) for re_idx in idx_tau_list ]

plt.plot(x,y, 'coral',linewidth=0.5)

for i,j,z in zip(x_label, y_label, label):
plt.text( i, j+0.03, z, ha='center', va='center', fontsize=8, color='steelblue')

pop_a = mpatches.Patch(color='coral', label=r'Cloud Effective Radius ($\mu m$)')
pop_b = mpatches.Patch(color='steelblue', label='Cloud Optical Thickness')

plt.legend(handles=[pop_a,pop_b],loc=2,fontsize=8)

if phase == 'liquid': plt.title("MODIS Liquid Cloud LUT",fontsize=8)
if phase == 'ice': plt.title("MODIS Ice Cloud LUT",fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[x_band_idx]) + r'\mu m$'
plt.xlabel(label,fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[y_band_idx]) + r'\mu m$'
plt.ylabel(label,fontsize=8)

plt.xlim(0,1.2)
plt.ylim(0,0.8)

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

plt.show()
plt.close()

### Plot liquid and ice LUT together

x_band_idx = 0
y_band_idx = 4

fig = plt.figure()

x = [reflectanceWater_data[index_re,x_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]
y = [reflectanceWater_data[index_re,y_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]

plt.plot(x,y,'--', color='coral',linewidth=0.5)

x = [reflectanceIce_data[index_re,x_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]
y = [reflectanceIce_data[index_re,y_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]

plt.plot(x,y,'--', color= 'steelblue',linewidth=0.5)

x = [reflectanceWater_data[:,x_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]
y = [reflectanceWater_data[:,y_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]

plt.plot(x,y,'--', color= 'coral',linewidth=0.5)

x = [reflectanceIce_data[:,x_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]
y = [reflectanceIce_data[:,y_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]

plt.plot(x,y,'--', color= 'steelblue',linewidth=0.5)

pop_a = mpatches.Patch(color='coral', label=r'Liquid Cloud')
pop_b = mpatches.Patch(color='steelblue', label='Ice Cloud')

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

plt.title("MODIS Cloud LUT",fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[x_band_idx]) + r'\mu m$'
plt.xlabel(label,fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[y_band_idx]) + r'\mu m$'
plt.ylabel(label,fontsize=8)

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

plt.show()
plt.close()

### Source Code

from pyhdf.SD import SD, SDC
from scipy import interpolate

import matplotlib.pyplot as plt
import numpy as np
import math
import matplotlib.patches as mpatches
import matplotlib.cm as cm

print( file.info() )

#----------------------------------------------------------------------------------------#
# Get SDS info

datasets_dic = file.datasets()

for idx,sds in enumerate(datasets_dic.keys()):
print( idx,sds )
sds_obj = file.select(sds)
print( sds_obj.info() )
print()
print()

#----------------------------------------------------------------------------------------#
# Get the Data

reflectanceWater = file.select('reflectanceWater')
reflectanceIce = file.select('reflectanceIce')
WaveLengths = file.select('WaveLengths')
OpticalThickness = file.select('OpticalThickness')

reflectanceWater_data = reflectanceWater.get()
reflectanceIce_data = reflectanceIce.get()
WaveLengths_data = WaveLengths.get()
OpticalThickness_data = OpticalThickness.get()

liq_lut_shape = reflectanceWater_data.shape
ice_lut_shape = reflectanceIce_data.shape

#----------------------------------------------------------------------------------------#
# Plot Nakajima and King LUT:

# band idx:
# 0 --> 0.65899998
# 1 --> 0.86500001
# 2 --> 1.24000001
# 3 --> 1.63999999
# 4 --> 2.13000011
# 5 --> 3.75

x_band_idx = 0
y_band_idx = 4

phase = 'ice' # liquid or ice

fig = plt.figure()

if phase == 'liquid':
x = [reflectanceWater_data[index_re,x_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]
y = [reflectanceWater_data[index_re,y_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]

if phase == 'ice':
x = [reflectanceIce_data[index_re,x_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]
y = [reflectanceIce_data[index_re,y_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]

plt.plot(x,y, 'steelblue',linewidth=0.5)

if phase == 'liquid':
x = [reflectanceWater_data[:,x_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]
y = [reflectanceWater_data[:,y_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]

if phase == 'ice':
x = [reflectanceIce_data[:,x_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]
y = [reflectanceIce_data[:,y_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]

plt.plot(x,y, 'coral',linewidth=0.5)

pop_a = mpatches.Patch(color='coral', label=r'Cloud Effective Radius ($\mu m$)')
pop_b = mpatches.Patch(color='steelblue', label='Cloud Optical Thickness')

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

if phase == 'liquid': plt.title("MODIS Liquid Cloud LUT",fontsize=8)
if phase == 'ice': plt.title("MODIS Ice Cloud LUT",fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[x_band_idx]) + r'\mu m$'
plt.xlabel(label,fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[y_band_idx]) + r'\mu m$'
plt.ylabel(label,fontsize=8)

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

plt.show()
plt.close()

#----------------------------------------------------------------------------------------#
# Add Re and Tau labels

#print(reflectanceIce_data)
#print(OpticalThickness_data)

x_band_idx = 0
y_band_idx = 4

phase = 'ice' # liquid or ice

fig = plt.figure()

if phase == 'liquid':
x = [reflectanceWater_data[index_re,x_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]
y = [reflectanceWater_data[index_re,y_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]

x_label_re = [reflectanceWater_data[index_re,x_band_idx,liq_lut_shape[2]-1] for index_re in np.arange(liq_lut_shape[0]) ]
y_label_re = [reflectanceWater_data[index_re,y_band_idx,liq_lut_shape[2]-1] for index_re in np.arange(liq_lut_shape[0]) ]

if phase == 'ice':
x = [reflectanceIce_data[index_re,x_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]
y = [reflectanceIce_data[index_re,y_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]

x_label_re = [reflectanceIce_data[index_re,x_band_idx,ice_lut_shape[2]-1] for index_re in np.arange(ice_lut_shape[0]) ]
y_label_re = [reflectanceIce_data[index_re,y_band_idx,ice_lut_shape[2]-1] for index_re in np.arange(ice_lut_shape[0]) ]

plt.plot(x,y, 'steelblue',linewidth=0.5)

for i,j,z in zip(x_label_re, y_label_re,label_re):
plt.text( i+0.04, j, z, ha='center', va='center', fontsize=8, color='coral')

if phase == 'liquid':
x = [reflectanceWater_data[:,x_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]
y = [reflectanceWater_data[:,y_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]

idx_tau_list = [0,3,6,9,12,15,18,21,34]

x_label = [reflectanceWater_data[0,x_band_idx,index_tau] for index_tau in idx_tau_list ]
y_label = [reflectanceWater_data[0,y_band_idx,index_tau] for index_tau in idx_tau_list ]

label = [str(OpticalThickness_data[re_idx]) for re_idx in idx_tau_list ]

if phase == 'ice':
x = [reflectanceIce_data[:,x_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]
y = [reflectanceIce_data[:,y_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]

idx_tau_list = [0,3,6,9,12,15,18,21,34]

x_label = [reflectanceIce_data[0,x_band_idx,index_tau] for index_tau in idx_tau_list ]
y_label = [reflectanceIce_data[0,y_band_idx,index_tau] for index_tau in idx_tau_list ]

label = [str(OpticalThickness_data[re_idx]) for re_idx in idx_tau_list ]

plt.plot(x,y, 'coral',linewidth=0.5)

for i,j,z in zip(x_label, y_label, label):
plt.text( i, j+0.03, z, ha='center', va='center', fontsize=8, color='steelblue')

pop_a = mpatches.Patch(color='coral', label=r'Cloud Effective Radius ($\mu m$)')
pop_b = mpatches.Patch(color='steelblue', label='Cloud Optical Thickness')

plt.legend(handles=[pop_a,pop_b],loc=2,fontsize=8)

if phase == 'liquid': plt.title("MODIS Liquid Cloud LUT",fontsize=8)
if phase == 'ice': plt.title("MODIS Ice Cloud LUT",fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[x_band_idx]) + r'\mu m$'
plt.xlabel(label,fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[y_band_idx]) + r'\mu m$'
plt.ylabel(label,fontsize=8)

plt.xlim(0,1.2)
plt.ylim(0,0.8)

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

plt.show()
plt.close()

#----------------------------------------------------------------------------------------#
# Plot liquid and ice LUT together

x_band_idx = 0
y_band_idx = 4

fig = plt.figure()

x = [reflectanceWater_data[index_re,x_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]
y = [reflectanceWater_data[index_re,y_band_idx,:] for index_re in np.arange(liq_lut_shape[0])]

plt.plot(x,y,'--', color='coral',linewidth=0.5)

x = [reflectanceIce_data[index_re,x_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]
y = [reflectanceIce_data[index_re,y_band_idx,:] for index_re in np.arange(ice_lut_shape[0])]

plt.plot(x,y,'--', color= 'steelblue',linewidth=0.5)

x = [reflectanceWater_data[:,x_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]
y = [reflectanceWater_data[:,y_band_idx,index_tau] for index_tau in np.arange(liq_lut_shape[2])]

plt.plot(x,y,'--', color= 'coral',linewidth=0.5)

x = [reflectanceIce_data[:,x_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]
y = [reflectanceIce_data[:,y_band_idx,index_tau] for index_tau in np.arange(ice_lut_shape[2])]

plt.plot(x,y,'--', color= 'steelblue',linewidth=0.5)

pop_a = mpatches.Patch(color='coral', label=r'Liquid Cloud')
pop_b = mpatches.Patch(color='steelblue', label='Ice Cloud')

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

plt.title("MODIS Cloud LUT",fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[x_band_idx]) + r'\mu m$'
plt.xlabel(label,fontsize=8)

label = r'MODIS Band $\lambda = ' + str(WaveLengths_data[y_band_idx]) + r'\mu m$'
plt.ylabel(label,fontsize=8)

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

plt.show()
plt.close()

Licence

Activity