Create MODIS C6 L2 Cloud Effective Radius histograms with python

There are several different ways to create an histogram of MODIS C6 L2 Cloud Effective Radius and filter the data.

Read SDS data and import needed modules

To download a MODIS granule from NASA lads server using python and ftp, go here. To learn how to read an HDF file with python go here

#!/usr/bin/env python

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 = './inputs/myd06_granules_viz/MYD06_L2.A2008015.1435.006.2013342100940.hdf'

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

file = SD(file_name, SDC.READ)

cer_id = file.select('Cloud_Effective_Radius')
cer_pcl_id = file.select('Cloud_Effective_Radius_PCL')
qa_id = file.select('Quality_Assurance_1km')
cpop_id = file.select('Cloud_Phase_Optical_Properties')
cm_id = file.select('Cloud_Mask_1km')

cer_data = cer_id.get()
cer_pcl_data = cer_pcl_id.get()
qa_data = qa_id.get()
cpop_data = cpop_id.get()
cm_data = cm_id.get()

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

cloud_mask_flag = bits_stripping(1,2,cm_data[:,:,0])
clear_sky_restoral_flag = bits_stripping(6,2,qa_data[:,:,3])

Plot Cloud Effective Radius histogram

Create a simple histogram with no filtering (meaning all pixels: cloudy or not; successful retrievals or not). It is just a baseline histogram:

attributes_cer = cer_id.attributes()

data = cer_data * attributes_cer['scale_factor']
#data = np.ravel(data)

data[ data < 0.0 ] = 0.0

data2 = data[ data > 0.0 ]
data3 = data[ data == 0.0 ]

plt.hist(data2, bins=50,color="lightblue", label='retrievals')
plt.hist(data3, bins=1,color="red", alpha=0.5, label='no retrievals')

plt.title('Cloud Effective Radius Histogram')
plt.xlabel(r'Cloud Effective Radius Retrievals at $2.1\mu m$')
plt.grid()
plt.legend(loc='upper right')
plt.xlim(-1.0,60.0)

plt.savefig("myd06_c6_cer_21_histogram.png", 
            bbox_inches='tight', dpi=200)
plt.show()
plt.close()

Plot Cloud Effective Radius histogram (successful retrievals only)

To keep only pixels where Cloud Effective Radius retrievals are successful, we just need to add a simple filter: data1[data1 > 0.0]

data = cer_data * attributes_cer['scale_factor']
data = np.ravel(data)
data = data[data > 0.0]

plt.hist(data, bins=50,color="lightblue")

plt.title('Cloud Effective Radius Histogram \n (successful retrievals only)')
plt.xlabel(r'Cloud Effective Radius Retrievals at $2.1\mu m$')
plt.grid()

plt.savefig("myd06_c6_cer_21_successful_histogram.png", 
            bbox_inches='tight', dpi=200)
plt.show()
plt.close()

Plot Cloud Effective Radius histogram (PCL; successful retrievals only)

attributes_cer_pcl = cer_pcl_id.attributes()

data = cer_pcl_data * attributes_cer_pcl['scale_factor']
data = np.ravel(data)
data = data[data > 0.0]

plt.hist(data, bins=50,color="lightblue")

plt.title('Cloud Effective Radius Histogram \n (PCL; successful retrievals only)')
plt.xlabel(r'Cloud Effective Radius Retrievals at $2.1\mu m$')
plt.grid()

plt.savefig("myd06_c6_cer_21_pcl_successful_histogram.png", 
            bbox_inches='tight', dpi=200)
plt.show()
plt.close()

Plot Cloud Effective Radius histogram (PCL and non-PCL; successful retrievals only)

data1 = cer_data * attributes_cer['scale_factor']
data1 = np.ravel(data1)
data1 = data1[data1 > 0.0]

data2 = cer_pcl_data * attributes_cer_pcl['scale_factor']
data2 = np.ravel(data2)
data2 = data2[data2 > 0.0]

data = np.concatenate((data1, data2), axis=0)

plt.hist(data, bins=50,color="lightblue")

plt.title('Cloud Effective Radius Histogram \n (PCL + non-PCL; successful retrievals only)')
plt.xlabel(r'Cloud Effective Radius Retrievals at $2.1\mu m$')
plt.grid()

plt.savefig("myd06_c6_cer_21_pcl_and_non_pcl_successful_histogram.png", 
            bbox_inches='tight', dpi=200)
plt.show()
plt.close()

same overlapping

data1 = cer_data * attributes_cer['scale_factor']
data1 = np.ravel(data1)
data1 = data1[data1 > 0.0]

data2 = cer_pcl_data * attributes_cer_pcl['scale_factor']
data2 = np.ravel(data2)
data2 = data2[data2 > 0.0]

data = np.concatenate((data1, data2), axis=0)

plt.hist(data1, bins=50,color="lightblue", label='non-PCL')
plt.hist(data2, bins=50,color="red", alpha=0.5, label='PCL')

plt.title('Cloud Effective Radius Histogram \n (PCL + non-PCL; successful retrievals only)')
plt.xlabel(r'Cloud Effective Radius Retrievals at $2.1\mu m$')
plt.grid()
plt.legend(loc='upper right')

plt.savefig("myd06_c6_cer_21_pcl_and_non_pcl_successful_overlapping_histogram.png", 
            bbox_inches='tight', dpi=200)
plt.show()
plt.close()

Plot Cloud Effective Radius histogram (PCL and non-PCL; successful retrievals only, cloud phase)

ice

data1 = cer_data * attributes_cer['scale_factor']
data1 = data1[ (data1 > 0.0) & (cpop_data == 3)]

data2 = cer_pcl_data * attributes_cer_pcl['scale_factor']
data2 = data2[ (data2 > 0.0) & (cpop_data == 3)]

data = np.concatenate((data1, data2), axis=0)

plt.hist(data, bins=50,color="lightblue")

plt.title('Ice Cloud Effective Radius Histogram \n (PCL + non-PCL; successful retrievals only)')
plt.xlabel(r'Cloud Effective Radius Retrievals at $2.1\mu m$')
plt.grid()

plt.savefig("myd06_c6_ice_cer_21_pcl_and_non_pcl_successful_histogram.png", 
            bbox_inches='tight', dpi=200)
plt.show()
plt.close()

liquid

data1 = cer_data * attributes_cer['scale_factor']
data1 = data1[ (data1 > 0.0) & (cpop_data == 2)]

data2 = cer_pcl_data * attributes_cer_pcl['scale_factor']
data2 = data2[ (data2 > 0.0) & (cpop_data == 2)]

data = np.concatenate((data1, data2), axis=0)

plt.hist(data, bins=50,color="lightblue")

plt.title('Liquid Cloud Effective Radius Histogram \n (PCL + non-PCL; successful retrievals only)')
plt.xlabel(r'Cloud Effective Radius Retrievals at $2.1\mu m$')
plt.grid()

plt.savefig("myd06_c6_liquid_cer_21_pcl_and_non_pcl_successful_histogram.png", 
            bbox_inches='tight', dpi=200)
plt.show()
plt.close()

both

data1 = cer_data * attributes_cer['scale_factor']
data1 = data1[ (data1 > 0.0) & (cpop_data == 3)]

data2 = cer_pcl_data * attributes_cer_pcl['scale_factor']
data2 = data2[ (data2 > 0.0) & (cpop_data == 3)]

data = np.concatenate((data1, data2), axis=0)

plt.hist(data, bins=50,color="lightblue", label='ice')

data1 = cer_data * attributes_cer['scale_factor']
data1 = data1[ (data1 > 0.0) & (cpop_data == 2)]

data2 = cer_pcl_data * attributes_cer_pcl['scale_factor']
data2 = data2[ (data2 > 0.0) & (cpop_data == 2)]

data = np.concatenate((data1, data2), axis=0)

plt.hist(data, bins=25,color="red", alpha=0.5, label='liquid')

plt.title('Cloud Effective Radius Histogram \n (PCL + non-PCL; successful retrievals only; cloud phase)')
plt.xlabel(r'Cloud Effective Radius Retrievals at $2.1\mu m$')
plt.grid()
plt.legend(loc='upper right')

plt.savefig("myd06_c6_cer_21_pcl_and_non_pcl_successful_cloud_phase_histogram.png", 
            bbox_inches='tight', dpi=200)
plt.show()
plt.close()

Plot Cloud Effective Radius histogram (using Cloud Mask Flag and CSR flag)

data1 = cer_data * attributes_cer['scale_factor']
data1 = data1[ (cloud_mask_flag < 2 ) & (clear_sky_restoral_flag == 0)]

data1[ data1 < 0.0 ] = 0.0

data2 = data1[ data1 > 0.0 ]
data3 = data1[ data1 == 0.0 ]

plt.hist(data2, bins=50,color="lightblue", label='successful')
plt.hist(data3, bins=1,color="red", alpha=0.5, label='not successful')

plt.title('Liquid Cloud Effective Radius Histogram \n (Cloud Mask = Cloudy, CSR = not restored)')
plt.xlabel(r'Cloud Effective Radius Retrievals at $2.1\mu m$')
plt.grid()
plt.legend(loc='upper right')
plt.xlim(-1.0,60.0)

plt.savefig("myd06_c6_liquid_cer_21_successful_and_not_histogram.png", 
            bbox_inches='tight', dpi=200)
plt.show()
plt.close()

References