How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?

Published: December 08, 2024
Updated: December 08, 2024
| Tags: CloudSat;
| Comments
DMCA.com Protection Status

Published: December 08, 2024

Updated: December 08, 2024

Tags: CloudSat;

DMCA.com Protection Status

Introduction

This guide provides a comprehensive, step-by-step tutorial for visualizing the CloudSat 2B-CLDCLASS-LIDAR Product using Python. The visualization includes important cloud attributes such as CloudLayerBase, CloudLayerTop, and CloudPhase, enabling an insightful representation of cloud structure and thermodynamic phase.

By the end of this tutorial, you'll be able to generate plots for cloud layers and phases from the CloudSat dataset, leveraging Python's powerful data processing and visualization libraries.

Table of contents

Prerequisites

Before proceeding, ensure you have the following Python libraries installed. You can install them using pip if they are not already available:

  • pyhdf: For reading HDF files (Hierarchical Data Format)
  • matplotlib: For data visualization
  • numpy: For numerical computations and data manipulation

Loading the Data

The CloudSat 2B-CLDCLASS-LIDAR product is stored in HDF format. You will load this data and extract relevant datasets for further analysis.

1
2
3
4
5
6
7
8
9
from pyhdf.SD import SD, SDC 
from pyhdf.HDF import *
from pyhdf.VS import *

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

Define the path to the dataset

1
2
3
4
5
path_to = '/Users/bmarchant/\
Desktop/MyDrive/Datasets/\
Courses/Introduction_to_python_for_satellite_remote_sensing/'

file_name = path_to + '2008183012329_11573_CS_2B-CLDCLASS-LIDAR_GRANULE_P_R04_E02.hdf'

Latitude and Longitude Extraction

Extract latitude and longitude information, which is crucial for geolocating cloud attributes.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
f = HDF(file_name, SDC.READ) 
vs = f.vstart()

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

latitude_table = Latitude[:]

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

Opening the HDF File

Load the file and inspect its contents to identify the datasets it contains.

1
2
3
4
5
file = SD(file_name, SDC.READ)

file_info = file.info()

print(file_info)  # Displays the number of datasets and metadata entries

Output

1
(13, 10)

Getting CloudLayerBase Product

CloudLayerBase represents the base altitude of cloud layers in kilometers. Extract and inspect its attributes.

1
2
3
4
5
6
7
sds_obj = file.select('CloudLayerBase')  # Select the dataset
CloudLayerBase = sds_obj.get()          # Retrieve data
sds_info = sds_obj.info()               # Get dataset info

print(sds_info)  # Dataset name, rank, dimensions, etc.
print(CloudLayerBase.shape)  # Shape of the data array
pprint.pprint(sds_obj.attributes())  # Print metadata attributes

Attributes Example:

1
2
3
4
5
6
7
{'_FillValue': -99.0,
'factor': 1.0,
'missing': -99.0,
'missop': '==',
'offset': 0.0,
'units': 'km',
'valid_range': [0.0, 20.0]}

CloudLayerBase first 10 rows example:

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

Output

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
[6.790016, 11.950026, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[6.7300158, 11.950026, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[6.790016, 12.010026, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[6.2500124, 11.950026, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[6.2500124, 12.010026, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[6.310013, 12.010026, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[6.2500124, 12.010026, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[6.2500124, 12.010026, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[2.4699984, 6.0100107, 12.010026, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[2.4699984, 6.2500124, 12.010026, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]

Getting CloudLayerTop Product

Similar to CloudLayerBase, CloudLayerTop provides the top altitude of cloud layers.

1
2
3
4
5
6
7
sds_obj = file.select('CloudLayerTop') # select sds

CloudLayerTop = sds_obj.get()

sds_info = sds_obj.info()

pprint.pprint( sds_obj.attributes() )

Attributes Example:

1
2
3
4
5
6
7
{'_FillValue': -99.0,
'factor': 1.0,
'missing': -99.0,
'missop': '==',
'offset': 0.0,
'units': 'km',
'valid_range': [0.0, 20.0]}

CloudLayerTop first 10 rows example:

1
2
for i in range(10):
    print([CloudLayerTop[i,j] for j in range(10)])

Output

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
[7.030018, 15.670052, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[7.0900183, 15.670052, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[7.0900183, 15.670052, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[7.1500187, 15.730052, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[7.1500187, 15.730052, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[7.1500187, 15.730052, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[7.1500187, 15.730052, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[7.1500187, 15.730052, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[2.6499982, 7.1500187, 15.730052, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]
[2.6499982, 7.1500187, 15.730052, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0, -99.0]

Getting CloudPhase Product

CloudPhase indicates the thermodynamic phase of the clouds (e.g., ice, liquid, or mixed phase).

1
2
3
4
5
6
7
sds_obj = file.select('CloudPhase') # select sds

CloudPhase = sds_obj.get()

sds_info = sds_obj.info()

pprint.pprint( sds_obj.attributes() )

Attributes Example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
{'_FillValue': -9,
'factor': 1.0,
'missing': -9,
'missop': '==',
'offset': 0.0,
'valid_range': [0, 2]}


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

Print CloudPhase first 10 rows example:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
[3, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[3, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[3, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[2, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[2, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[2, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[2, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[2, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[3, 2, 1, 0, 0, 0, 0, 0, 0, 0]
[3, 2, 1, 0, 0, 0, 0, 0, 0, 0]

Plotting Cloud Layers

Visualize the cloud layers, using their base and top altitudes.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
cldclass_lidar_start_idx = 000
cldclass_lidar_end_idx = 1000

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

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

ax.set_facecolor('xkcd:lightblue')

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

plt.legend()

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

xticks_idx = [int(i) for i in np.linspace(cldclass_lidar_start_idx,cldclass_lidar_end_idx,10)]
xticks_latitude = [ "{0:.2f}".format(latitude_table[i][0]) for i in xticks_idx]

plt.xticks(xticks_idx, xticks_latitude,rotation=45)

plt.xlabel('Latitude')
plt.ylabel('Altitude (in km)')

plt.title('Cloud Mask')

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

plt.close()

How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?
How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?

Plotting Cloud Phase

Distinguish cloud phases using different colors for visualization.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
f = plt.figure()
ax = f.add_subplot(111)

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

ax.set_facecolor('xkcd:lightblue')

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

plt.legend()

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

xticks_idx = [int(i) for i in np.linspace(cldclass_lidar_start_idx,cldclass_lidar_end_idx,10)]
xticks_latitude = [ "{0:.2f}".format(latitude_table[i][0]) for i in xticks_idx]

plt.xticks(xticks_idx, xticks_latitude,rotation=45)

plt.xlabel('Latitude')
plt.ylabel('Altitude (in km)')

plt.title('Cloud thermodynamic phase')

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

plt.close()

How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?
How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?

Vectorization

The code above performs well and is sufficient for plotting a small portion of a CloudSat track. However, if you want to speed up the process, you can optimize the code by vectorizing it. Although vectorized code can be more complex to understand and relies on Python features like broadcasting, it significantly improves performance.

To vectorize our code, we can create a simple grid or heatmap and populate it with data. Below is an example that demonstrates how to plot the first cloud layer, regardless of the cloud phase:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
y_bin_size = 0.1
y_max = 20.0

y_bin_count = int(y_max/y_bin_size)

CloudLayerBase_Indexes =  ( CloudLayerBase / y_bin_size ).astype(int)

CloudLayerTop_Indexes =  ( CloudLayerTop / y_bin_size ).astype(int)


x = np.arange(0, y_bin_count)
y = np.arange(0, 37081)

xx, yy = np.meshgrid(x, y)


heatmap = ((CloudLayerBase_Indexes[:,0].reshape(37081,1)  <= xx) 
        & (xx <= CloudLayerTop_Indexes[:,0].reshape(37081,1)  ))



plt.imshow(heatmap[0:1000,:].T,aspect='auto',origin='lower')

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

How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?
How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?

How to plot the first ice cloud layer:

1
2
3
4
5
6
7
8
heatmap = ( ( CloudLayerBase_Indexes[:,0].reshape(37081,1)  <= xx )
        & ( xx <= CloudLayerTop_Indexes[:,0].reshape(37081,1) )
        & ( np.repeat( (CloudPhase[:,0].reshape(37081,1) == 1)  ,y_bin_count,axis=1)  ) )

plt.imshow(heatmap[0:1000,:].T,aspect='auto',origin='lower')

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

How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?
How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?

Now, let's iterate over each layer and cloud phase.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
y_bin_size = 0.1
y_max = 20.0

y_bin_count = int(y_max/y_bin_size)

CloudLayerBase_Indexes =  ( CloudLayerBase / y_bin_size ).astype(int)

CloudLayerTop_Indexes =  ( CloudLayerTop / y_bin_size ).astype(int)

heatmap = np.zeros((37081, y_bin_count))

x = np.arange(0, y_bin_count)
y = np.arange(0, 37081)

xx, yy = np.meshgrid(x, y)

for phase in [1,2,3]:
    for layer in range(10):
        heatmap_sub = ( ( CloudLayerBase_Indexes[:,layer].reshape(37081,1)  <= xx )
                & ( xx <= CloudLayerTop_Indexes[:,layer].reshape(37081,1) )
                & ( np.repeat( (CloudPhase[:,layer].reshape(37081,1) == phase),y_bin_count,axis=1)  ) ).astype(int)

        heatmap_sub[ heatmap_sub == 1 ] = phase

        heatmap = heatmap + heatmap_sub

Plotting cloud phase

1
2
3
4
plt.imshow(heatmap[0:1000,:].T,aspect='auto',origin='lower',interpolation='nearest')

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

How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?
How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?

Visualization Enhancements:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
cmap =  [(0.6784313725490196, 0.8470588235294118, 0.9019607843137255, 1.0)] + \
         [(1.0, 1.0, 1.0, 1.0)] + \
         [(0.9803921568627451, 0.5019607843137255, 0.4470588235294118, 1.0)] + \
         [(0.0, 0.0, 1.0, 1.0)]


cmap = mpl.colors.ListedColormap(cmap)

bounds = [0.0, 1.0, 2.0, 3.0, 4.0]

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

img = plt.imshow(heatmap[0:1000,:].T,
                 aspect='auto',
                 origin='lower',
                 interpolation='nearest',
                 cmap=cmap, 
                 norm=norm)

cbar_bounds = [0.0, 1.0, 2.0, 3.0, 4.0]
cbar_ticks = [0.5,1.5,2.5,3.5]  
cbar_labels = ['Clear', 'Ice Cloud','Mixed','Liquid Cloud']

cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_labels, fontsize=10)


plt.xlim(cldclass_lidar_start_idx,cldclass_lidar_end_idx)

xticks_idx = [int(i) for i in np.linspace(cldclass_lidar_start_idx,cldclass_lidar_end_idx,10)]
xticks_latitude = [ "{0:.2f}".format(latitude_table[i][0]) for i in xticks_idx]

plt.xticks(xticks_idx, xticks_latitude,rotation=45)



yticks_idx = [int(i/y_bin_size) for i in np.linspace(0,y_max,20)]
yticks = [int(i) for i in np.linspace(0,y_max,20)]

plt.yticks(yticks_idx, yticks )


plt.xlabel('Latitude')
plt.ylabel('Altitude (in km)')

plt.title('Cloud thermodynamic phase')



plt.grid(linestyle='--')

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

How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?
How to Plot CloudSat 2B-CLDCLASS-LIDAR Product Using Python ?

References

Links Site
CloudSat 2B-CLDCLASS-LIDAR Product ProcessDescription and Interface Control Document cloudsat.cira.colostate.edu
Creative Commons License

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

Author profile-image

Benjamin

Hi, I'm Ben, the Founder of moonbooks.org. I work as a research scientist specializing in Earth satellite remote sensing, particularly focusing on Fire detection using VIIRS and ABI sensors. My areas of expertise include Python, Machine Learning, and Open Science. Feel free to reach out if you notice any errors on the site. I am constantly striving to enhance it. Thank you!

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

of