• The objectives of this notebook is to provide examples about working with CloudSat L2 and AUX files, including:
    • Background information of CALIPSO Project, its functionality and mission;
    • How to access to the CALIPSO datasets via NASA's **Atmospheric Science Data Center**;
    • How to draw a satellite orbit;
    • How to find and segment the orbit.
  • The files required for this notebook including:
  • Files above are CALIPSO L2 Lidar Cloud Profiles and Layers Product, Granule Number 02702, observed on Oct 30 2006, Which consist with what we discussed in 01_MODIS_L1B.ipynb.
  • *.hdf is the HDF4 file, Section 2 introduced the way to get these files.
  • *.h5 is the HDF5 file which can be converted from HDF4 files. See 01_MODIS_L1B.ipynb and Section 3 for more information.

In [0]:
__author__ = 'ATSC-301 UBC'


About CALIPSO Project

**Wikipedia *CALIPSO* Page**

CALIPSO is a joint NASA (USA) and CNES (France) environmental satellite, built in the Cannes Mandelieu Space Center, which was launched atop a Delta II rocket on April 28, 2006. Its name stands for Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations.

Passive and active remote sensing Instruments on board the CALIPSO satellite monitor aerosols and clouds 24 hours a day. CALIPSO is part of the "A Train", flying in formation with several other satellites (Aqua, Aura and CloudSat).

You can findout more about CALIPSO's news and applications via **NASA's Satellite Mission Page**.

Access to the CALIPSO datasets

CALIPSO datasets are available in NASA's **Atmospheric Science Data Center**.

In order to get the datasets, we do the following steps:

  • We choose Cloud Profiles and Layers Product in CALIPSO L2 Lidar product family, it contains the following parameters:

    • Integrated Attenuated Backscatter
    • Column Reflectance
    • Lidar Depolarization Ratio
    • Cloud Base Height
    • Cloud Top Height
    • Cloud Midlayer Temperature
  • There are three different choice of cloud layer resolution, we select **1/3km version**.

  • Click Order Data and search available granules based on geolocation and time:
  • Click "Search for Granules" and add your ideal files to the "Shopping Cart":
  • Click "View" and send your order:
  • On order page, you are required to provide personal informations. If you are previously regeistered in any NASA's Earthdata system, you can also use that account to send your order. Order processed very quick and you will receive an email about what to do next.

Convert HDF4 files to HDF5

After you downloaded the CALIPSO files via FTP client. You need to convert HDF4 files to HDF5.

We have discussed this part in 01_MODIS_L1B.ipynb. But here we have a large number of files, a recommended way is to use h4toh5convert in batch sctipt.

  • For Windows user: a *.bat batch script runs in DOS Terminal is located in _share/. copy *.bat to the direcctory of your HDF4 files (here say _data/CloudSat_L2), open the Terminal, cd to that path and type:
  • For Linux user: there is a *.sh script locates in the same place, similarly, you do the steps above and type:
bash BASH_batch_h4toh5.sh
  • I do not have a Mac but I think Bash works for Mac <--------------- !!!

External tools


Except for what we will do in the following sections, **ccplot** is another choice for you:

ccplot is an open source command-line program for plotting profile, layer and earth view data sets from CloudSat, CALIPSO and Aqua MODIS products.

ccplot works for Windows, Linux and Mac. Visualising Data from CloudSat and CALIPSO Satellites - Peter Kuma.pdf in \_docs is its documentation. It also has a **Github Page**.


There is a CALIPSO_tools provided by Paul Cottle in _lib/, it organzied based on Python pandas module and can findout the cloud top height in Vertical Feature Mask, VFM product in CALIPSO L2 Lidar product family. Currenly, the problem it is a bit slow, and actually we do not use VFM product in this notebook, so we provide the tool, but we will not use it.

If you are interested in this tool, here is the example codes for how to use it:

# import the module
import sys
sys.path.insert(0, '_libs/CALIPSO_tools')
# search the file
Cal_obj = Calipso.Calipso(h5file = hdf5_VFM[0], raw = True)
# Segment the orbit
lats = [45, 50]; lons = []; alts = [0,40000]
Cal_obj.window_select(latrange = lats,lonrange = lons,altrange = alts, inplace = True)
# Calculate Cloud top
Cal_masked = Cal_obj.feature_mask('cloud')
Cal_cloudtops = Calipso.findtops(Cal_masked)
# A simple visulation
plt.plot(Cal_cloudtops.VFM.index.get_level_values('Lats').values, Cal_cloudtops.VFM.values)

Import modules and tools

In [0]:
import h5py
import glob
import math
import scipy.io
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
% matplotlib inline

In order to read the time of orbit, here we need Python module datetime

In [0]:
import datetime

Read CALIPSO through h5py

Searching all *.h5 files in _data/CALIPSO_L2/ directory, Make sure that you have the required HDF5 files exist in it.

In [0]:
print("CALIPSO 333m Cloud Layer files found:\n{}".format(hdf5_calipso))

Get HDF5 file object

In [0]:
calipso_obj=h5py.File(hdf5_calipso[0], 'r')

Read different parameters

We have introduced how to play with HDF5 in 01_MODIS_L1B.ipynb, so here we skip that part.

Geolocation field

In [0]:

We need to make sure that longitude is monotonous decrease, because satellite is possible to cross the International Date Line (180E/180W line)

In [0]:
for id in range(0, len(lon)-1):
    if lon[id+1] > lon[id]:
        lon[id+1] = lon[id+1]-360

Currently NOT sure about height in the cell below, Lidar_data_Altitude in metadata has a shape of 583L, but the size of datasets are (55926L, 5L). I guess Lidar_data_Altitude represents the top height of every 55920/583 = 95 data points, but I don't how to make it corresponds to 5L in vertical dimension. <----------- !!!

In [0]:
height=calipso_obj['metadata'][0][-5] # <----------------- Not Sure, use the code below instead
height=np.array([0, 1, 2, 3, 4])

Lidar Cloud top Height

In [0]:
cloud_top_height[cloud_top_height == -9999]=np.nan

Orbit time

In [0]:

View the orbit track

Now we try to plot the orbit track of the CALIPSO file based on lon, lat and profile_time

We need to know when the orbit start, and it is the 2nd element of metadata (You can check it via .keys() or other tools)

In [0]:

Then we convert the date string to datetime object.

In [0]:
orbitStart=datetime.datetime.strptime(tai_start_string, '%Y-%m-%dT%H:%M:%S.%fZ')

Create another datetime object which includes all time info. of data points.

In [0]:
#python datetime objects in utc
for the_time in profile_time:
    date_time=orbitStart + datetime.timedelta(seconds=float(the_time))

Finally, we plot the orbit with time on a map. There are more than 50000 data points in the file, so we have to skip some of them.

The code is a bit complex, here is the key idea:

time_step=1500; # plot every 1500 data points
label_step=2; # label every 2 ploted data points
for i in np.arange(0, len(date_day), time_step):
    # every timestep
    proj.plot(x[i], y[i], 'bo', markersize=8)  # plot the data points
    # labelling 
    if(i % (time_step*label_step) == 0):  # just like "mod" in MATLAB 
        time_string=date_day[i].strftime('%H:%M UCT')   # .strftime() is a call of output format for datetime object
        axis.text(x_str[i], y_str[i], '{0:s}'.format(time_string), \     # axis.text(x, y, 'string', **keywords)
                  fontsize=8, fontweight='bold', ha='left', \
                  bbox={'edgecolor':'w', 'facecolor':'w', 'alpha':0.875, 'pad':0})

In [0]:
time_step=1500; label_step=2;
Orbit_number=calipso_obj['metadata'][0][10] # <---- Granule Number when orbit start
# locations of time label 
lat_str=lat-1; lon_str=lon+4 # lat/lon for labels
# figure
fig=plt.figure(figsize=(15, 15))
proj = Basemap(projection='mill', resolution='c', \
    llcrnrlat=-90, urcrnrlat=90, \
    llcrnrlon=-240, urcrnrlon=120)
proj.drawcoastlines(color=[0.375, 0.375, 0.375])
# draw parallels and meridians.
proj.drawparallels(np.arange(-60, 90, 30), labels=[1, 0, 0, 0])
proj.drawmeridians(np.arange(0, 360, 60), labels=[0, 0, 0, 1])
# compute native x,y coordinates of grid.
x, y=proj(lon, lat)
x_str, y_str=proj(lon_str, lat_str)
# line of track
proj.plot(x, y, linewidth=3.5, color='k', linestyle='-')
# =================================== main part ================================= #
for i in np.arange(0, len(date_day), time_step):
    # every timestep
    proj.plot(x[i], y[i], 'bo', markersize=8)
    # labelling 
    if(i % (time_step*label_step) == 0):
        time_string=date_day[i].strftime('%H:%M UCT')
        axis.text(x_str[i], y_str[i], '{0:s}'.format(time_string), \
                  fontsize=8, fontweight='bold', ha='left', \
                  bbox={'edgecolor':'w', 'facecolor':'w', 'alpha':0.875, 'pad':0})
# ================================================================================ #
# starting point in red
proj.plot(x[0], y[0], 'ro', markersize=8)
# bold axis
[i.set_linewidth(2) for i in axis.spines.itervalues()]
# title
title_str='CALIPSO Track\nGranule Number: 0' + str(Orbit_number)
axis.set_title(title_str, fontsize=14, fontweight='bold')
# savefig
plt.savefig('_figures/02_CALIPSO_track.png', dpi=450, facecolor='w', edgecolor='w',
            orientation='portrait', papertype='a4', format='png',
            transparent=True, bbox_inches='tight', pad_inches=0,

Segment the orbit track

In [0]:

Read the data we saved in 01_MODIS_L1B.ipynb.

Remembering that we use

scipy.io.savemat('_share/01_MODIS_L1B_TBright', {'longitude': longitude, 'latitude': latitude, 'cloud_Tbright': cloud_Tbright})

to save the file.

In [0]:

Mask the NaN value in Tbring

In [0]:
cloud_Tbright_masked=np.ma.masked_where(np.isnan(cloud_Tbright), cloud_Tbright)

Plot the track with clouds

We plot the Brightness Temperature to represents the clouds, and plot the orbit track. The code is similar to the last part of 01_MODIS_L1B.ipynb.

The orbit is a line, but here I plot it with super big linewidth to make it clear.

In [0]:
lonlim=[-155, -127]
latlim=[30, 60]
# Colormap
CMap=plt.cm.gray_r#gist_heat # hot, afmhot, gnuplot
proj=Basemap(resolution='l', projection='lcc', \
            lat_1=30, lat_2=60, lat_0=45, lon_0=-140, \
            llcrnrlon=-155, llcrnrlat=30, \
            urcrnrlon=-110, urcrnrlat=56)
# create figure, add axes
fig=plt.figure(figsize=(12, 12))
## parallels and meridians.
parallels=np.arange(-90, 90, 5)
meridians=np.arange(0, 360, 5)
proj.drawparallels(parallels, labels=[1, 0, 0, 0],\
                  fontsize=10, latmax=90)
proj.drawmeridians(meridians, labels=[0, 0, 0, 1],\
                  fontsize=10, latmax=90)
# draw coast 
proj.drawcoastlines(linewidth=1.5, linestyle='solid', color=[0.25, 0.25, 0.25])
# compute native x,y coordinates of grid.
x, y=proj(lon_L1B, lat_L1B)
x_track, y_track=proj(lon, lat)
# pcolor plot
CS=proj.pcolor(x, y, cloud_Tbright_masked, cmap=CMap, vmin=210, vmax=275)
# plot the track
proj.plot(x_track, y_track, linewidth=15, color='r', linestyle='-', alpha=0.625)
# title
ax.set_title('Cloud and Orbit Track',\
             fontweight='bold', fontsize=14)
# Save figure
plt.savefig('_figures/02_Cloud_and_Orbit_Track.png', dpi=450, facecolor='w', edgecolor='w',
            orientation='portrait', papertype='a4', format='png',
            transparent=True, bbox_inches='tight', pad_inches=0,
# Show

Define segmented region based on latitude

Based on the figure above, we can segment the orbit by latitude:

$35^\circ N < Latitude < 52^\circ N$

We use `numpy.searchsorted` to segment. This function has been discussed in 01_MODIS_L1B.ipynb.

Here, latitude is the "bins", and latitude=35, latitude=50 is the value.

In [0]:
minlat_indices=np.searchsorted(lat.flat, 35, 'right')-1
maxlat_indices=np.searchsorted(lat.flat, 52, 'right')-1
# Cut

Since the longitude of the orbit is monotone decreasing according to our discussion.

But the output of np.searchsorted is naturally in a monotone increasing sense, so here we change the "direction" of lon_seg.

In [0]:
lon_seg=lon_seg[::-1, :];
#cloud_top_height=cloud_top_height[::-1, :];

We save the result as *.mat

In [0]:
time_str=[0] * len(date_day_seg)
for i in range(len(date_day_seg)):
    time_str[i]=date_day_seg[i].strftime('%H:%M UCT')

scipy.io.savemat('_share/02_CALIPSO_L2', {'longitude': lon_seg, 'latitude': lat_seg, 'time': time_str, 'surf_elev': surf_elev_seg, 'Cloud_Top_Height': cloud_top_height})

Plot the result

We plot the CALIPSO lidar cloud top height on the segmented orbit. Essentially, it is:

plt.plot(lat_seg, cloud_top_height)

Mask the NaN value.

In [0]:
cloud_top_height_masked=np.ma.masked_where(np.isnan(cloud_top_height), cloud_top_height)

In [0]:
from mpl_toolkits.axes_grid1 import host_subplot
import mpl_toolkits.axisartist as AA
# figure
fig=plt.figure(figsize=(10, 5))
# host_subplot
host = host_subplot(111, axes_class=AA.Axes)
# two more axis
par2 = host.twiny()
par3 = host.twiny()
# Still not so sure about how this works, currently looks good : )
offset = -25
new_fixed_axis = par2.get_grid_helper().new_fixed_axis
par2.axis["top"] = new_fixed_axis(loc="bottom", axes=par2, offset=(0, offset))

offset = -62.5
new_fixed_axis = par3.get_grid_helper().new_fixed_axis
par3.axis["top"] = new_fixed_axis(loc="bottom", axes=par3, offset=(0, offset))
# axis label <---- don't know why fontsize=14, fontweight='bold' doesn't work
host.set_ylabel('Lidar Cloud top Height (km)', fontsize=14, fontweight='bold')
host.set_xlabel('\n\nGeolocation (longitude/latitude)')
par3.set_xlabel('Orbit time')
# axis lim
par2.set_xlim(lat_seg.min(), lat_seg.max())
par3.set_xlim(0, 1)
host.set_xlim(lon_seg.min(), lon_seg.max())
host.set_ylim(-0.5, 10)
# plot
# fill the surface elevation
host.fill_between(lon_seg.flat[:], surf_elev_seg.flat[:], baseline, \
                  where=surf_elev_seg.flat[:]>=baseline, facecolor=[0.5, 0.5, 0.5], interpolate=False)
host.plot(lon_seg, cloud_top_height_masked[:, 0], linewidth=1.5, linestyle='-', color='k')
# xtick and ytick
host.set_yticks([1, 3, 5, 7, 9])
host.set_xticks(np.linspace(lon_seg.min(), lon_seg.max(), 6))
par2.set_xticks(np.linspace(lat_seg.min(), lat_seg.max(), 6))
count=0; labels_str=[0] * 6
for time_id in np.linspace(0, len(date_day_seg)-1, 6).astype(int):
    labels_str[count]=date_day_seg[time_id].strftime('%H:%M UCT')
host.set_title('CALIPSO Lidar Cloud top Height', fontsize=16, fontweight='bold')
# savefig
plt.savefig('_figures/02_CALIPSO_Lidar_Cloud_top_Height.png', dpi=450, facecolor='w', edgecolor='w',
            orientation='portrait', papertype='a4', format='png',
            transparent=True, bbox_inches='tight', pad_inches=0,

In [0]:

In [0]: