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

Published: June 27, 2019
DMCA.com Protection Status

Published: June 27, 2019

DMCA.com Protection Status

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

Table of contents

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

file = SD('./MODIS_C6_LUTS/examples/ocean_PY_C6_reflectance_sol_08_view_03_relazm_13_wspeed_06pt00_2017151142552.hdf', SDC.READ)

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')
ParticleRadiusWater = file.select('ParticleRadiusWater')
reflectanceIce = file.select('reflectanceIce')
ParticleRadiusIce = file.select('ParticleRadiusIce')
WaveLengths = file.select('WaveLengths')
OpticalThickness = file.select('OpticalThickness')

reflectanceWater_data = reflectanceWater.get()
ParticleRadiusWater_data = ParticleRadiusWater.get()
reflectanceIce_data = reflectanceIce.get()
ParticleRadiusIce_data = ParticleRadiusIce.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

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

# 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

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

#print(ParticleRadiusWater_data)
#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]) ]

    label_re = [str(ParticleRadiusWater_data[re_idx]) for re_idx in np.arange(ParticleRadiusWater_data.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]) ]

    label_re = [str(ParticleRadiusIce_data[re_idx]) for re_idx in np.arange(ParticleRadiusIce_data.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

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

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

file = SD('./MODIS_C6_LUTS/examples/ocean_PY_C6_reflectance_sol_08_view_03_relazm_13_wspeed_06pt00_2017151142552.hdf', SDC.READ)

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')
ParticleRadiusWater = file.select('ParticleRadiusWater')
reflectanceIce = file.select('reflectanceIce')
ParticleRadiusIce = file.select('ParticleRadiusIce')
WaveLengths = file.select('WaveLengths')
OpticalThickness = file.select('OpticalThickness')

reflectanceWater_data = reflectanceWater.get()
ParticleRadiusWater_data = ParticleRadiusWater.get()
reflectanceIce_data = reflectanceIce.get()
ParticleRadiusIce_data = ParticleRadiusIce.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(ParticleRadiusWater_data)
#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]) ]

    label_re = [str(ParticleRadiusWater_data[re_idx]) for re_idx in np.arange(ParticleRadiusWater_data.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]) ]

    label_re = [str(ParticleRadiusIce_data[re_idx]) for re_idx in np.arange(ParticleRadiusIce_data.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()

References

Links Site
Determination of the Optical Thickness and Effective Particle Radius of Clouds from Reflected Solar Radiation Measurements. Part I: Theory .ametsoc.org
Overview of the MODIS Collection 6 Cloud Optical Property (MOD06) Retrieval Look-up Tables nasa.gov
Look Up Tables (LUTs) nasa.gov
Creative Commons License

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Author profile-image

Benjamin

Developer in my spare time, I’m interested in Earth observation and satellite-based fire detection. Passionate about Python, machine learning, and open science.

Skills
  • Remote sensing
  • Python
  • Machine Learning
  • Tensorflow 2.0
Image

of